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Summary 


There  has  been  increasing  interest  in  developing  numerical  methods  and  accurate 
physical  models  for  solving  computational  fluid  dynamics  (CFD)  problems  of  hypersonic 
continuum  and  rarefied  flows.  In  this  work,  we  establish  a  framework  for  validating  nu¬ 
merical  methods  and  physical  models  employed  in  popular  CFD  codes.  The  first  main 
objective  of  this  work  is  to  assess  the  ability  of  current  state-of-the-art  methods  to  simu¬ 
late  challenging  hypersonic  flow  problems  by  comparing  to  well  characterized  experiments. 
The  second  main  objective  is  to  provide  benchmark  numerical  solutions  to  hypersonic 
double-cone  laminar  flows  that  can  be  used  as  code  validation  cases.  This  report  has  been 
derived  from  Dr.  Ioannis  Nompelis’  Ph.  D.  Thesis  at  the  University  of  Minnesota. 

We  have  simulated  experiments  of  hypersonic  laminar  double-cone  flows  that  were 
performed  at  the  Large  Energy  National  Shock  (LENS)  facility.  These  experiments  were 
conducted  as  part  of  a  validation  effort  for  continuum  and  particle-based  codes.  The 
double-cone  flow  was  chosen  because  it  exhibits  strong  viscous/inviscid  and  shock  interac¬ 
tions,  and  earlier  work  had  shown  that  it  is  very  challenging  to  compute  and  it  is  sensitive 
to  real-gas  effects.  The  initial  set  of  experiments  consisted  of  low  enthalpy  hypersonic 
flows,  and  our  simulations  of  these  flows  were  aiming  to  evaluate  the  quality  of  the  nu¬ 
merics.  Subsequently,  following  completion  of  the  work  on  low  enthalpy  experiments,  new 
double-cone  flow  experiments  at  high  enthalpy  were  performed.  The  high  enthalpy  exper¬ 
iments  will  be  used  for  validating  and  improving  the  physical  models  and  nonequilibrium 
chemical  reaction  rates.  Both  the  low  enthalpy  and  high  enthalpy  LENS  experiments 
were  designed  to  result  in  laminar  flows,  such  that  uncertainty  in  the  simulations  due  to 
turbulence  modeling  is  not  present. 

The  initial  comparisons  between  numerical  predictions  and  experimental  measure¬ 
ments  for  low  enthalpy  experiments  exhibited  discrepancies.  The  simulations  predicted 
correctly  the  size  of  the  region  where  the  flow  was  separated  as  well  as  the  surface  pres¬ 
sure.  However,  the  heat  transfer  rate  to  the  first  cone  was  over-predicted.  We  investigated 
the  source  of  this  discrepancy  by  simulating  the  entire  experiment,  including  the  flow  in 
the  facility  hypersonic  nozzle.  We  found  that  vibrational  energy  freezes  at  the  throat 
conditions  in  the  free-stream  for  these  low  density  nitrogen  experiments.  This,  and  weak 
vibrational  energy  accommodation  at  the  surface  due  to  surface  slip  effects  was  responsible 
for  the  discrepancies  observed.  We  also  considered  the  effect  of  weak  flow  non-uniformity  in 
the  test-section  and  finally  obtained  complete  agreement  with  the  experimental  measure¬ 
ments.  This  is  an  important  step  towards  validation  of  numerical  methods  for  two  reasons. 
First,  we  have  obtained  agreement  to  well  within  the  experimental  uncertainty  for  the  most 
carefully  designed  and  most  highly  instrumented  experiments  known  to  us.  Second,  we 
have  developed  a  methodology  of  characterizing  the  LENS  experiments  in  terms  of  the 


state  of  the  free-stream  using  CFD. 

We  investigated  the  effects  of  numerics  on  the  double-cone  flow,  and  in  particular 
we  focused  on  the  amount  of  numerical  dissipation  exhibited  by  the  different  numerical 
methods  employed.  We  found  that  the  predicted  size  of  the  separation  zone  on  a  given  mesh 
is  directly  related  to  the  amount  of  numerical  dissipation  in  the  numerical  flux  evaluation 
method.  We  also  conducted  grid  convergence  studies  for  different  numerical  schemes  and 
showed  that  by  making  comparisons  to  our  computed  results,  a  direct  assessment  of  the 
grid-convergence  rate  of  a  given  numerical  method  can  be  made.  It  is  important  to  note 
that  a  complete  grid  convergence  study  is  not  required  for  this  assessment  to  be  made. 

We  also  simulated  high  enthalpy  double-cone  experiments  of  nitrogen  and  air.  We 
obtained  good  agreement  for  nitrogen,  particularly  at  low  free-stream  Reynolds  number. 
Agreement  with  the  air  flows  was  only  qualitative.  We  performed  nozzle  flow  simulations 
but  we  were  unable  to  improve  agreement  for  the  air  flow,  and  this  is  probably  due  to 
uncertainty  in  the  reaction  rates.  This  is  an  area  where  further  research  should  be  done. 

Finally,  we  conducted  a  study  of  experiments  performed  using  the  newly  developed 
LENS-X  expansion  tunnel.  Experiments  at  high  enthalpy  air  performed  using  the  expan¬ 
sion  tunnel  will  result  in  free-stream  with  realistic  composition,  and  thus  will  help  improve 
nonequilibrium  air  chemistry  models. 


Introduction 


1.  Motivations 

There  has  been  a  great  deal  of  interest  in  developing  prototype  re-entry  vehicles  and 
hypersonic  air-breathing  vehicles  for  access  to  space.  These  vehicles  will  be  designed  to 
operate  in  the  high  enthalpy  regime,  where  the  flow  exhibits  high  enthalpy  or  “real-gas”  ef¬ 
fects.  Real-gas  effects  are  triggered  by  the  high  temperatures  encountered  in  high  enthalpy 
flows.  These  effects  include  excitation  of  the  internal  energy  modes  of  the  gas,  as  well  as 
dissociation  of  polyatomic  species  and  ionization.  These  physical  phenomena  occur  over  a 
range  of  time  scales,  and  the  flow  is  frequently  in  chemical  and  thermal  nonequilibrium. 
Additionally,  depending  on  the  shape  of  the  vehicle  and  the  flight  conditions,  shock/shock 
interactions  and  viscous-inviscid  interactions  will  likely  occur  during  flight  at  hypersonic 
speeds.  The  combined  effects  of  these  phenomena  and  real-gas  chemistry  may  dramatically 
alter  the  flight  characteristics  of  the  vehicle.  Therefore,  all  of  these  effects  must  be  taken 
into  consideration  in  the  development  stage  of  future  hypersonic  vehicles. 

The  development  and  design  of  such  vehicles  is  aided  by  the  use  of  experimentation 
and  numerical  simulation.  There  are  difficulties  associated  with  both  approaches.  It  is 
difficult  to  produce  ground-based  hypersonic  flow  experiments  that  are  realistic  in  terms 
of  the  enthalpy  and  composition  of  the  gas.  It  is  also  difficult  to  perform  accurate  nu¬ 
merical  simulations,  as  the  physics  involved  is  not  fully  understood.  For  example,  in  the 
sudden  expansion  of  high  temperature  air,  the  coupled  vibrational  relaxation  and  chemical 
recombination  process  is  poorly  understood.  This  is  an  example  where  modeling  of  the 
thermochemistry  is  important,  and  uncertainty  is  introduced  into  the  simulation  through 
the  chemical  reaction  rates.  It  also  appears  that  critically  important  hypersonic  flow 
phenomena  such  as  shock/shock  and  viscous-inviscid  interactions  cannot  be  predicted  ac¬ 
curately.  Hypersonic  flows  with  regions  of  separation  induced  by  shock/shock  interactions 
have  been  very  difficult  problems  to  tackle  computationally.  These  flow  phenomena  are 
associated  with  a  large  range  of  length  scales,  particularly  at  high  Reynolds  numbers.  In 
many  cases,  a  large  number  of  grid  points  is  required  in  order  to  accurately  predict  the 
basic  flow  structure.  This  increases  substantially  the  cost  of  numerical  simulations. 

It  is  clear  that  although  there  are  different  aspects  of  hypersonic  flow  -  each  exhibiting 
its  own  idiosyncrasies  -  the  physical  processes  may  not  be  separated  and  decoupled  in  nu¬ 
merical  simulations  for  practical  problems.  Therefore  it  is  important  to  get  a  perspective 
view  on  the  reliability  and  versatility  of  numerical  prediction  techniques  used  in  real-life 
applications.  The  accuracy  of  popular  numerical  methods  used  in  computational  fluid  dy¬ 
namics  (CFD)  codes  must  be  assessed,  along  with  the  accuracy  of  existing  nonequilibrium 
chemistry  models.  Furthermore,  it  is  important  to  understand  the  limitations  of  the  partic- 


1 


ular  numerical  methods  and  physical  models.  To  achieve  this,  we  typically  require  that  the 
numerical  simulation  is  tested  against  a  problem  for  which  the  solution  is  known.  For  the 
complex  problems  of  interest,  analytical  solutions  cannot  be  obtained.  As  an  alternative, 
numerical  prediction  techniques  are  required  to  reproduce  experimental  measurements  at 
well  characterized  conditions. 

2.  Double-Cone  Flows 

Flows  generated  by  double-cone  geometries  placed  in  hypersonic  free-streams  have 
been  studied  in  the  past.  Earlier  studies  of  these  flows  focused  primarily  on  testing  and 
improving  nonequilibrium  chemistry  models.  This  is  because  the  double-cone  flows  are 
sensitive  to  high  enthalpy  effects,  and  it  has  been  speculated  that  their  sensitivity  to  real- 
gas  chemistry  would  allow  for  discriminating  between  air  chemistry  models.  The  double¬ 
cone  flows  are  also  very  difficult  to  simulate  numerically,  and  thus  they  can  be  utilized  as 
a  comprehensive  test  of  CFD  codes. 

The  flow  generated  by  a  double-cone  configuration  placed  in  a  hypersonic  free-stream, 
exhibits  strong  viscous-inviscid  and  shock/shock  interactions.  The  oblique  shock  generated 
by  the  first,  shallow-angle  cone,  interacts  with  the  shock  generated  by  the  second  cone. 
When  the  second  cone  angle  is  larger  than  a  critical  angle,  this  shock  is  detached,  and  the 
interaction  is  very  strong.  A  schematic  of  the  flow  generated  by  a  double-cone  in  a  nitrogen 
free-stream  at  a  Mach  number  of  about  12  is  shown  in  Fig.  1.  In  this  figure,  the  model 
has  a  first  cone  half- angle  of  25°  and  a  second  cone  half-angle  of  55°.  The  shock  structure 
and  subsonic  regions  are  outlined,  and  the  shear  layers  are  visualized.  The  oblique  shock 
interacts  with  the  bow  shock  generated  by  the  second  cone,  and  a  transmitted  shock  is 
formed,  which  impinges  on  the  surface  downstream  of  the  cone-cone  juncture.  The  adverse 
pressure  gradient  at  the  cone-cone  juncture  causes  the  boundary  layer  to  separate.  The 
separated  region  that  forms  at  the  corner  creates  its  own  shock,  which  interacts  with  the 
bow  shock.  This  shifts  the  point  of  interaction,  which  in  turn  alters  the  separation  zone. 
This  process  feeds  back  on  itself  until  the  flow  reaches  steady  state  -  if  a  steady  state 
exists.  A  supersonic  jet  is  formed  along  the  surface  of  the  second  cone  downstream  of  the 
impingment  point,  and  it  undergoes  a  series  of  isentropic  compressions  and  expansions. 
The  flow  behind  the  strong  bow  shock  is  subsonic,  and  therefore  the  shape  of  the  jet 
influences  the  shock  structure. 

The  complex  flow  phenomena  exhibited  by  double-cone  geometries  under  hypersonic 
conditions  make  the  flow  very  difficult  to  compute  with  CFD.  This  is  mainly  because  of 
two  reasons.  First,  there  is  a  large  range  of  scales  that  must  be  resolved.  In  particular,  the 
boundary  layer  that  grows  along  the  first  cone  influences  the  location  of  the  oblique  shock. 
A  large  number  of  grid  points  with  stretching  near  the  wall  is  required  to  correctly  resolve 
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Figure  1.  Schematic  diagram  of  a  typical  flow-field  generated  by  a  25°-55°  double¬ 
cone  geometry  placed  in  a  Mach  12  free-stream.  The  flow  is  from  left  to  right,  and 
a  slice  is  shown  of  flow  that  is  circularly  symmetric  about  the  axis  of  the  geometry. 


the  boundary  layer  growth  near  the  tip,  the  scale  of  which  is  much  smaller  compared  to 
the  size  of  the  geometry.  Also,  the  triple  point  and  the  point  of  re-attachment  must  be 
correctly  resolved,  because  they  directly  influence  the  size  of  the  separation  zone.  Thus,  fine 
grid  spacing  is  also  required  near  the  point  of  interaction  and  re-attachment.  Therefore,  a 
very  fine  mesh  is  required  for  double-cone  simulations.  Second,  the  presence  of  separation 
results  in  a  feedback  mechanism  by  which  the  separation  zone  continuously  evolves,  and  this 
process  takes  a  long  time  to  reach  steady  state  relative  to  the  characteristic  flow  time.  We 
refer  to  a  characteristic  flow  time  as  the  time  it  takes  for  a  particle  to  travel  the  length  of  the 
geometry  if  it  is  moving  at  the  free-stream  velocity.  Therefore,  the  combined  requirements 
of  having  high  grid  resolution  and  having  to  simulate  large  amounts  of  physical  time,  make 
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such  simulations  extremely  expensive.  We  may  seek  to  reduce  the  number  of  grid  points  by 
making  use  of  higher  order  low  dissipation  methods.  Even  in  this  case,  an  efficient  implicit 
time-integration  method  is  required  in  order  to  have  reasonable  turn-around  times  on  a 
parallel  computer.  Therefore,  despite  the  simplicity  of  the  geometry,  the  double-cone  flow 
is  a  stringent  test  for  numerical  methods  and  algorithms. 

Under  high  enthalpy  conditions,  the  chemical  processes  that  take  place  in  the  flow 
are  tightly  coupled  with  the  fluid  motion  through  energy  exchange  mechanisms.  In  high 
enthalpy  double-cone  flows,  dissociation  takes  place  both  in  the  high  temperature  region 
behind  the  bow  shock,  and  inside  the  separation  zone.  In  these  regions,  the  highly  ener¬ 
getic  particles  undergo  collisions  and  lose  some  of  their  energy  by  causing  dissociation  of 
polyatomic  species.  By  this  process,  energy  is  removed  from  the  bulk  flow  and  stored  in 
the  form  of  chemical  energy.  The  rate  at  which  this  process  takes  place  greatly  impacts 
the  mean  flow.  Thus,  the  size  of  the  embedded  region  of  separation  is  sensitive  to  chem¬ 
ical  effects  in  the  flow.  Additionally,  the  size  of  the  separation  zone  may  be  measured 
experimentally  by  surface  measurements  such  as  surface  pressure  and  heat  transfer  rate. 
Thus,  we  can  make  direct  comparisons  between  numerical  predictions  and  experimental 
measurements.  This  is  a  desirable  attribute  for  a  flow  that  may  be  used  as  a  test  case 
for  validating  and  improving  nonequilibrium  chemistry  models.  Therefore,  provided  that 
a  given  numerical  method  is  able  to  accurately  predict  double-cone  flows  in  the  absence  of 
flow-field  chemistry,  thermochemical  models  can  be  tested  using  double-cone  flows  of  high 
enthalpy. 

3.  Review  of  Related  Work 

Hypersonic  double-cone  and  double-wedge  flows  have  been  studied  extensively  both 
numerically  and  experimentally.  In  the  1990’s,  attempts  were  made  to  compare  simulations 
results  with  experimental  data  for  these  interacting  flows.  The  results  of  these  studies 
showed  good  agreement  in  some  cases,  but  in  many  cases  the  efforts  were  not  as  successful. 

Olejniczak  et  al.  (1999)  made  use  of  double-wedges  placed  in  a  hypersonic  free-stream 
at  the  California  Institute  of  Technology  Free-Piston  Shock  Tunnel  T5  to  test  different 
nonequilibrium  chemistry  models.  They  did  not  obtain  good  agreement  in  terms  of  the 
separation  zone  size  and  distribution  of  surface  quantities,  and  they  attributed  this  to 
inadequate  modeling  of  the  thermochemistry  and  flow  physics.  Around  the  same  time, 
Olejniczak  et  al.  (1997a)  studied  inviscid  shock  interactions  on  double-wedges  and  found 
that  in  many  cases  the  gas  dynamics  have  a  dominant  effect  on  the  surface  pressure  dis¬ 
tribution.  This  explained  the  pressure  variations  near  the  re-attachment  point  of  the  flow, 
which  were  thought  to  be  due  to  real-gas  effects.  Attempts  to  reproduce  double-cone  exper¬ 
iments  in  high  enthalpy  nitrogen  were  also  unsuccessful  [Olejniczak  et  al.  (1997b)].  These 


4 


efforts  were  aimed  at  testing  nonequilibrium  chemistry  models  using  hypersonic  flows  with 
shock  interactions. 

More  recent  simulations  of  double-cone  flows  by  Wright  et  al.  (2000)  of  experiments 
performed  at  the  Princeton  University  Mach  8  Wind  Tunnel  showed  good  agreement  with 
measurements.  They  showed  that  laminar  simulations  predicted  the  surface  quantities  at 
low  Reynolds  numbers.  Also,  for  higher  Reynolds  numbers,  where  the  flow  was  believed 
to  be  transitional,  turbulent  simulations  showed  good  agreement. 

4.  Scope  of  the  Current  Study 

In  this  work,  we  use  computational  fluid  dynamics  to  study  hypersonic  double-cone 
flows.  Double-cone  flows  under  both  low  enthalpy  and  high  enthalpy  conditions  are  inves¬ 
tigated.  This  work  aims  toward  the  validation  of  physical  models  and  numerical  methods 
that  are  employed  in  computational  fluid  dynamics  codes.  It  serves  as  a  starting  point  for 
code  validation,  and  is  intended  to  understand  the  extent  to  which  numerical  methods  can 
predict  hypersonic  low-density  flows  with  laminar  regions  of  separation  in  the  presence  of 
chemical  reactions.  We  attempt  to  systematically  isolate  each  of  the  physical  aspects  that 
introduce  uncertainty.  For  example,  the  uncertainty  in  turbulence  modeling  is  eliminated 
by  ensuring  that  the  Reynolds  number  is  low  enough  such  that  the  flow  remains  laminar. 
Also,  high  enthalpy,  or  “real-gas”  effects  are  minimized  in  flows  at  low  enough  enthalpy, 
and  for  which  perfect-gas  numerical  predictions  can  be  made.  The  geometry  considered 
is  one  that  produces  strong  viscous-inviscid  and  shock/shock  interactions,  while  the  flow 
remains  steady.  We  compare  numerical  simulation  results  with  experimental  data  obtained 
in  the  Large  Energy  National  Shock  tunnel  (LENS)  facility,  at  the  Calspan  -  University  at 
Buffalo  Research  Center  (CUBRC). 

In  the  following  two  sections  of  this  report,  we  present  the  governing  equations  for  the 
continuum  flow  of  a  gas  in  vibrational  and  chemical  nonequilibrium.  The  Navier-Stokes 
equations  for  a  chemically  reacting  and  vibrationally  relaxing  air  mixture  are  presented,  as 
well  as  their  numerical  method  of  solution.  We  use  standard  numerical  techniques  which 
are  popular  among  researchers  and  have  been  widely  used  for  practical  applications.  In  the 
same  context,  the  standard  models  used  for  vibrational  relaxation  and  chemical  reaction 
rates  are  presented,  as  well  as  the  models  for  the  transport  coefficients. 

The  fourth  section  focuses  on  low  enthalpy  hypersonic  double-cone  flows  in  inert  envi¬ 
ronments.  Comparisons  between  experimental  measurements  and  blind  simulation  results 
are  presented,  with  the  aim  to  provide  an  objective  assessment  of  the  state  of  current 
simulation  methods.  The  fifth  section  is  concerned  with  the  numerical  aspects  of  these 
simulations,  and  presents  sensitivity  studies  related  to  the  discretization  of  the  equations 
and  the  accuracy  of  the  numerical  methods.  An  assessment  of  the  rate  of  convergence  is 
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also  made  for  this  particular  class  of  flows. 

The  sixth  section  contains  a  more  detailed  analysis  of  the  low  enthalpy  simulation 
results,  and  numerical  simulations  of  the  entire  experiments  are  presented.  These  include 
simulations  of  the  hypersonic  facility  nozzle,  where  comparisons  are  made  between  nu¬ 
merical  predictions  and  measurements  made  by  the  Pitot  pressure  probes  at  the  nozzle 
exit  plane.  As  a  result,  an  important  idiosyncrasy  of  these  low  density  flows  is  uncovered. 
Using  this  approach,  the  free-stream  conditions  are  revisited  and  additional  characteriza¬ 
tion  of  the  experimental  facility  is  done  using  CFD.  Further  double-cone  simulations  and 
comparisons  with  experimental  measurements  are  presented. 

The  seventh  section  discusses  work  on  high  enthalpy  double-cone  flows  that  was  per¬ 
formed  guided  by  the  methodology  developed  by  studying  the  low  enthalpy  experiments. 
Comparisons  between  experimental  measurements  and  numerical  predictions  for  high  en¬ 
thalpy  air  and  nitrogen  mixtures  are  presented. 

The  eighth  section  examines  high  enthalpy  experiments  performed  using  an  expansion 
tube  facility.  This  is  done  by  using  CFD  to  simulate  the  entire  experiment.  The  objective 
is  to  help  the  design  of  future  high  enthalpy  experiments  with  non-reacted  free-streams. 
Using  this  approach,  additional  characterization  of  the  facility  is  provided  in  terms  of 
available  test-time. 
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Mathematical  Formulation 

In  this  section,  we  introduce  the  set  of  coupled  partial  differential  equations  that  de¬ 
scribes  the  dynamics  of  the  flow-field.  For  the  high  enthalpy  flows  of  interest,  the  standard 
mass,  momentum  and  energy  equations  are  not  sufficient  to  describe  the  thermochemical 
state  of  the  gas  and  physics  involved.  The  high  temperatures  inherent  to  these  flows  trigger 
a  number  of  chemical  and  energy  relaxation  processes,  or  “real-gas”  effects,  and  the  flows 
are  frequently  in  substantial  nonequilibrium.  We  account  for  real-gas  effects  by  introducing 
finite-rate  chemical  reactions  and  internal  energy  relaxation  to  the  conservation  equations. 

The  Navier-Stokes  equation  set  has  been  extended  so  that  it  accurately  describes  the 
physics  of  a  gas  that  is  vibrationally  excited  and  chemically  reacting  [Candler  (1988), 
Olejniczak  (1997)],  and  it  is  only  briefly  presented  here.  The  equation  of  state  of  the  gas 
is  developed.  The  models  used  for  the  transfer  of  mass  between  different  chemical  species, 
transport  of  momentum,  and  the  transfer  of  energy  between  different  energy  modes  are 
discussed  along  with  their  underlying  assumptions. 

1.  Basic  Assumptions 

The  flowfields  are  assumed  to  be  accurately  described  by  a  continuum  formulation. 
For  this  assumption  to  be  valid,  the  Knudsen  number  defined  as  the  ratio  of  the  mean  free 
path  to  the  length  scale  of  the  flow  Kn  =  A/L,  must  be  small.  For  the  flows  of  interest  the 
Knudsen  number  is  of  the  order  of  0.1-0.001  based  on  the  model  size,  and  thus  the  flowfield 
is  in  the  continuum  regime.  It  is  implicit  in  the  continuum  formulation  that  there  is  little 
statistical  variation  at  any  point  in  the  flow  and  therefore  the  continuum  description  of 
the  viscous  fluxes  is  consistent. 

We  restrict  the  formulation  to  non-ionized,  chemically  reacting  flows.  The  rotational 
state  of  the  gas  typically  equilibrates  within  a  few  collisions  and  we  can  assume  that  the 
translational  and  rotational  states  of  the  gas  are  in  equilibrium  at  a  common  temperature. 
In  this  work  we  do  not  consider  electronic  excitation  of  the  particles. 

2.  The  Conservation  Equations 

For  a  chemically  reacting  mixture  of  ns  species,  of  which  nd  are  diatomic,  there  are 
ns  species  mass  conservation  equations  of  the  form 

+  V  •  (psu)  =  -V  •  ( psvs )  +  ws  , 

where  ps  is  the  density  of  species  s,  v$  is  the  species  diffusion  velocity,  and  ws  is  the  rate  of 
production  of  the  species  due  to  chemical  reactions.  By  summing  the  species  density  over 
all  species  we  get  the  mixture  density  p.  If  we  sum  the  species  mass  conservation  equation 
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over  all  species  s,  the  source  terms  that  appear  on  the  right  hand  side  sum  identically 
to  zero  and  we  recover  the  traditional  continuity  equation.  The  models  used  to  describe 
the  diffusion  velocities  and  the  form  of  the  source  terms  will  be  discussed  in  subsequent 
sections. 

The  conservation  of  linear  momentum  is  expressed  as 

d(pu)  ,  _  '  _  rX  _  _ 

^  -  +  V  •  ( pu  <g>  u  +  pi)  —  V  •  r  , 

where  p  is  the  thermodynamic  pressure  and  f  is  the  viscous  stress  tensor.  This  vector 
equation  represents  3  momentum  equations,  one  for  each  spatial  dimensions  of  the  prob¬ 
lem.  The  above  expression  assumes  that  there  are  no  external  body  forces  exerted  on 
the  particles.  Note  that  diffusion  velocities  do  not  appear  in  the  momentum  equation  as 
their  contributions  sum  identically  to  zero  (see  Olejniczak  (1997)  for  the  details  of  the 
calculation) . 

The  conservation  of  total  energy  of  the  mixture  is  given  by 

9E  ns^ 

~K7  +  V  •  (( E  +  p)u )  =  V  •  (f u)  —  V  •  (qt  4-  qr  +  qv)  —  V  •  ^(ps^s^s)  , 


where  (f ,  qr,  qv  are  the  translational,  rotational  and  vibrational  heat  flux  vectors,  and  hs 
is  the  total  enthalpy  per  unit  mass  of  species  s.  The  terms  present  in  the  energy  equation 
will  be  discussed  separately. 


In  our  formulation  we  assume  that  vibration- vibration  energy  coupling  is  very  strong 
such  that  all  vibrational  modes  are  in  equilibrium  with  each  other.  Thus  the  vibrational 
energy  of  the  diatomic  species  can  be  described  by  a  single  vibrational  temperature.  This 
is  a  significant  simplification  as  we  only  have  one  vibrational  energy  conservation  equation, 
while  in  principle  we  need  to  have  a  separate  equation  for  each  of  the  polyatomic  species. 
The  equation  for  the  vibrational  energy  per  unit  volume  is  given  as 


dEv 

dt 


nd 

+  V  •  ( Evu )  =  -V  •  qv  -  V  •  ^(psevsVs)  +  wv 

5=1 


5 


where  evs  is  the  vibrational  energy  per  unit  mass  of  species  s.  The  source  term  appearing 
on  the  right  hand  side  of  the  equation  will  be  discussed  in  the  subsequent  sections. 


3.  Equations  of  State 

In  this  section,  we  show  how  conserved  variables  relate  to  the  primitive  flow  variables, 
such  as  temperature,  and  pressure  that  appears  in  the  conservation  equations.  The  total 
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energy  per  unit  volume  of  the  gas  is  defined  as 

ns  j  nd  ns 

E  —  Y,PsCvsT  +  -p|«|2  +  evsPs  +  h°Ps  > 
s= 1  s=l  s=l 


where  Cvs  is  the  translational-rotational  specific  heat  at  constant  volume  for  species  s, 
and  is  given  by 

n  —  C  ( tr )  4-  C  ( rot ) 

Ws  —  Ws  i  Ws 

The  translational  and  rotational  specific  heats  are  given  by 


C  (tr)  = 


n  (.rot) 
uvs 


2  Ms 
R_ 
Ms 


s  <  nd  , 


where  R  is  the  universal  gas  constant  and  Ms  is  the  species  molecular  weight.  Thus  the 
first  term  in  the  energy  corresponds  to  the  internal  energy  of  the  gas  due  to  thermal 
translational-rotational  motion  of  particles.  The  second  term  represents  the  bulk  kinetic 
energy  of  the  gas.  The  third  term  represents  the  vibrational  energy  of  the  polyatomic 
species.  The  vibrational  energy  per  unit  mass  evs  of  species  s  is  related  to  the  vibrational 
temperature  Tv  assuming  a  Boltzmann  distribution  exists.  In  this  work,  we  use  a  simple 
harmonic  oscillator  model  to  represent  the  vibrational  excitation  of  diatomic  species  and 
this  internal  energy  is  expressed  as 


R  evs 

6vs  Ms  exp  (9VS/TV)  —  1  ’ 

where  0VS  is  the  characteristic  temperature  of  vibration  for  species  s.  Note  that  it  is  not 
possible  to  solve  this  equation  directly  for  Tv,  and  therefore  the  vibrational  temperature 
must  be  computed  iteratively  for  a  given  vibrational  energy  Ev  —  X^s=i  Psevs-  The  final 
term  in  the  expression  for  the  total  energy  represents  the  chemical  energy  of  formation 
associated  with  each  species. 


The  thermodynamic  pressure  p  is  related  to  the  translational-rotational  temperature 
through  a  perfect  gas  law  and  is  obtained  by  Dalton’s  Law  of  partial  pressures  of  the 
species  as 


R 


p= £?•= 


S=  1 


ns 

and  R  = 

$=1 


ps  R 

P  Ms  ’ 


The  enthalpy  per  unit  mass  of  species  s  is  defined  to  be 


4.  The  Diffusive  Terms 


In  this  section  we  describe  the  viscous,  thermal  and  mass  diffusion  fluxes  that  appear 
in  the  conservation  equations,  and  their  constitutive  relations. 

The  viscous  stress  tensor  is  formed  assuming  a  Newtonian  fluid  and  is  given  by 

f  =  /x(Vu  +  (Vu)T)  +  A(V  •  u)I  , 
with  the  Stokes  hypothesis  for  the  bulk  viscosity 

2/i  +  3A  —  0  , 


where  n  is  the  kinematic  viscosity  of  the  mixture.  The  viscosity  /i  depends  on  the  state  of 
the  fluid  and  its  composition  and  will  be  discussed  in  detail  below. 

The  heat  fluxes  are  assumed  to  obey  Fourier’s  Law  and  given  by 


Qi  —  , 


where  i  represents  any  of  the  internal  energy  modes  of  the  gas,  with  n *  and  T;  the  cor¬ 
responding  conductivity  and  temperature.  Since  the  translational  and  rotational  states 
of  the  gas  are  in  equilibrium  at  a  common  temperature,  it  is  convenient  to  combine  the 
corresponding  heat  fluxes  into  a  sum  of  the  conductivities  times  the  temperature  gradient 
as 

Q  —  Qt  +  Qr  =  ~{Kt  T  «r)VT  . 

The  thermal  conductivities  for  the  translational,  rotational,  and  vibrational  energy  modes 
are  determined  from  an  Eucken  relation  [Vincenti  and  Kruger  (1986)]  as 

=  \»sC^r)  ,  Krs  =  VsCv{rt)  ,  Kvs  =  »sC^ih)  , 

where  fis  and  ks  are  the  species  viscosities  and  thermal  conductivities.  The  vibrational 
specific  heat  is  computed  directly  using 


C1  (vib) 
Vs 


Mass  diffusion  is  driven  by  gradients  of  concentration,  pressure,  and  temperature. 
However,  for  the  flows  of  interest  only  the  terms  due  to  concentration  gradients  are  sig¬ 
nificant.  Thermal  diffusion  is  important  in  regions  where  gradients  are  large,  but  in  these 
regions  the  fundamental  continuum  assumption  is  typically  invalid  and  we  can  no  longer 
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apply  the  Navier-Stokes  equations  [Boyd  et  al.  (1995)].  The  baro-diffusion  terms  are  negli¬ 
gible.  Therefore,  for  this  work  the  mass  diffusion  is  solely  based  on  concentration  gradients 
and  Fick’s  Law  is  used.  The  mass  diffusion  flux  is 

psvs  =  —pDsV(—)  . 

P 

The  multicomponent  diffusion  coefficients  T>s  are  replaced  with  a  single  binary  diffusion 
coefficient  X>,  which  is  derived  assuming  a  constant  Lewis  number,  Le,  defined  as 


Le  =  X> 


pCp 


K 


where  Cp  is  the  translational-rotational  specific  heat  at  constant  pressure  and  k  is  the 
thermal  conductivity  of  the  mixture.  The  effects  of  multicomponent  diffusion  have  been 
neglected,  however  they  should  be  included  if  light  species,  such  as  hydrogen,  are  present 
in  the  flow. 


It  now  remains  to  provide  expressions  for  the  mixture  viscosity,  the  thermal  conduc¬ 
tivity  and  the  species  viscosities.  The  mixture  viscosity  p  and  thermal  conductivity  k  are 
obtained  from  the  species  viscosities  and  conductivities  ps,  ks  using  Wilke’s  semi-empirical 
mixing  rule  [Wilke  (1950)] 


p 


E 


XsPs 


«-E 


Xsf^>s 

4>s 


where 


<Ps  =  J2Xr 


1  + 


\Ps_ 

Pr 


Mr 

M, 


1/4.2  r 


-1 


and  cs  —  ps/ p  is  the  mass  fraction  of  species  s.  The  model  of  Blottner  is  used  for  the 
species  viscosities  ps,  which  are  calculated  from  curve  fits  of  the  form 


ps  =  0.1  exp[(yls  InT  +  Bs )  InT  +  Cs]  , 


where  As,  Bs,  and  Cs  are  constants.  This  model  is  valid  for  a  large  range  of  temperature 
(to  about  10,000  K),  and  is  suitable  for  this  work. 


5.  The  Chemical  Source  Terms 

Each  of  the  source  terms  ws  that  appear  in  the  mass  conservation  equations  represents 
the  rate  of  production  and  destruction  of  species  s  due  to  chemical  reactions.  In  this  section 
we  present  the  chemical  source  terms  that  appear  in  the  mass  conservation  equations. 
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In  the  high  enthalpy  flows  that  we  examine  in  this  work,  two  different  gases  are  used: 
nitrogen  and  air.  The  treatment  of  nitrogen  can  be  considered  as  a  special  case  of  the 
air  gas,  and  thus  the  set  of  conservation  equations  for  nitrogen  flows  is  a  subset  of  the 
equation  set  that  describes  air  flows.  For  high  temperature  air  with  no  ionization,  there 
are  five  chemical  reactions  that  are  important,  involving  five  species:  N2,  O2,  NO,  N,  and 
O.  The  chemical  reactions  involved  are 


N2  T  M  2N  +  M 
O2  T  M  20  -f-  M 
NO  +  M^N  +  O  +  M 
N2  +  O  NO  +  N 
NO  +  O  ^  02  +  N 


The  reactions  are  presented  so  that  they  are  endothermic  in  the  forward  direction.  The 
first  three  reactions  represent  dissociation  of  the  diatomic  species  N2,  02,  and  NO  due 
to  collisions.  In  these  reactions,  M  is  arbitrary  and  represents  the  collision  partner  which 
provides  the  energy  to  break  the  chemical  bond.  The  three  corresponding  reverse  reactions 
are  three-body  recombination,  during  which  the  collision  partner  M  absorbs  the  energy 
released  during  the  formation  of  chemical  bonds.  The  two  remaining  reactions  are  the 
Zeldovich  exchange  reactions,  which  are  the  primary  source  of  nitric  oxide  (NO)  in  a 
hypersonic  flow. 

The  rates  of  the  forward  and  backward  reactions  take  an  Arrhenius  form  governed  by 
the  forward  and  backward  rate  coefficients  kfrjn  and  /q>rm  respectively.  We  may  write  the 
rate  of  each  reaction  as  a  sum  of  the  forward  and  backward  rates 


Pn  Pn  Pm  _  t  Pn2  Pm 

MNMNMm~  flrnM„2  Mm. 

Po2  Pm 
"In.  Mm 


Til  =  K 

m 

'R'2  =  Y1  [kb2mM~'M~M  kf2 

m  L  lvlo  lvlo  lvlm  *** o2 

Pn  Po  Pm  u  Pno  Pm 


1Z3 = y:  h3 


m  Mn  Mm  "v°m  Mr, 


n  “o  *TJ-rn 

Pn2  Po 


‘no  ±v±m  - 


M»2Mq 


JKiNO 

t)  7„  Po-2.  Pn  Pno  Po 

~  kb5M~W  ~  kf5l \TW 

lvx02  J'JN  •‘“no  ■‘“o 
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The  chemical  source  terms  can  now  be  expressed  in  terms  of  the  individual  reaction  rates 


™n2 

=  MN2  (7£i  +  7^4) 

Wo2 

=  M02  (7Z2  —%s) 

™NO 

=  Mno  (ft3  —  7^-4  +  7^5) 

=  Mn  (—277i  —  77-3  —  774 

-775) 

wD 

=  M0  (—2772  —  773  +  774 

+  775) 

We  note  that  the  sum  of  all  the  source  terms  is  identically  zero,  so  that  conservation  of 
total  mass  in  satisfied,  and  elemental  conservation  of  atomic  nitrogen  and  oxygen  is  also 
satisfied. 

The  forward  and  backward  reaction  rate  coefficients  are  in  general  affected  by  the 
level  of  thermal  nonequilibrium  in  the  flow.  In  this  work  we  use  the  model  of  Park  (1986) 
for  vibration-dissociation  coupling.  The  Park  model  has  been  used  extensively,  since  it 
is  simple  to  implement  and  has  been  shown  to  give  satisfactory  results  for  a  variety  of 
flow  conditions  [Olejniczak  (1997)].  The  Park  model  assumes  that  the  reaction  rates  are 
functions  of  an  effective  temperature,  and  therefore  we  represent  the  forward  reaction  rate 
coefficients  by  the  modified  Arrhenius  form 

kfm  =  ^fmTei fT  exP(— ^m/3eff)  > 

where  C/  ,  r}m,  and  0m  are  found  from  curve  fits  to  experimental  data.  In  this  work,  the 
values  recommended  by  Park  (1988)  are  used.  The  effective  temperature  Teff  that  appears 
in  the  expression  is  some  function  of  the  translational-rotational  and  vibrational  temper¬ 
atures.  Park  recommends  that  the  geometric  average  of  the  translational  and  vibrational 
temperatures  is  used  for  the  forward  rates  of  the  dissociation  reactions,  as 

Te ff  =  VTTv  . 


The  backward  reaction  rate  coefficients  are  related  to  the  forward  reaction  rate  coefficients 
through  the  equilibrium  constant  defined  as  their  ratio 


K^JT)  = 


to  (~^eff  ) 

kbm 


The  equilibrium  constant,  Ke q,  is  also  a  curve  fit  to  experimental  data,  and  is  a  function 
only  of  T 

Keqm  Cm  exp(Aim  T  A.2m%  T  A.^mZ  A.^mZi  T  A^mZ  )  , 

where  Z  —  10,000/T.  The  constants  that  appear  in  the  expression  for  Ke<im  are  given  by 
Park  (1985). 
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6.  Vibrational  Energy  Source  Term 

The  source  term  that  appears  on  the  right  hand  side  of  the  vibrational  energy  equation 
(wv)  is  presented  here.  This  source  term  is  due  to  the  exchange  of  energy  between  the 
translational-rotational  and  vibrational  modes  of  the  gas.  Additionally,  the  creation  of 
diatomic  molecules,  due  to  recombination  reactions,  contributes  to  the  vibrational  energy 
source  term. 


There  are  many  energy  exchange  mechanisms  between  the  different  internal  energy 
modes  of  the  gas,  including  exchange  between  translational  and  vibrational  modes,  ex¬ 
change  between  rotational  and  vibrational  modes,  and  exchange  between  vibrational  modes 
of  the  separate  species.  In  this  work,  we  have  assumed  that  the  vibration-vibration  ex¬ 
changes  are  very  fast,  and  thus  the  vibrational  modes  of  the  gas  are  in  equilibrium  at  the 
same  vibrational  temperature  Tv.  The  vibration-translation  and  vibration-rotation  energy 
exchange  rates  are  combined  in  a  single  energy  exchange  rate  Qt-v,  and  the  Landau-Teller 
model  is  used  [Vincenti  and  Kruger  (1986)].  The  energy  exchange  rate  is  given  as 


ALT) 


Qt-v  s  ~  Ps~ 


•"Vs 


<TS> 


where  e*s  is  the  vibrational  energy  per  unit  mass  of  species  s  evaluated  at  the  translational 
temperature  (equilibrium),  and  <  rs  >  is  the  molar  averaged  Landau-Teller  relaxation 
time,  and  is  given  by  Lee  (1985)  as 


<  Ts  >- 


Er-^r 

Y.rXr/Tsr 


where  Xr  has  been  defined  previously  and  rsr  is  Landau-Teller  inter-species  relaxation 
time  given  by  Millikan  and  White  (1963),  as 


Tsr  =  -  expL4sr(T  1/3  —  0.015/i^4)  —  18.42]  ,  (p  in  atm) 
P 

Asr  =  1.16  x  KT3^2^3, 

Ms  Mr 

Msr  "  Ms  +  Mr  • 


Part  of  the  vibrational  energy  source  terms  is  due  to  the  creation  and  destruction  of 
diatomic  species  in  chemical  reactions.  For  this  work,  we  assume  that  the  molecules  that 
are  recombining  and  molecules  that  dissociate  carry  vibrational  energy  given  by  the  local 
vibrational  temperature.  Therefore,  the  rate  of  production  of  vibrational  energy  due  to 
chemical  reactions  is  just 

Qvs  ^s^vs  • 
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The  vibrational  energy  source  term  ws  is  simply  the  sum  of  the  contributions  of  the  two 
processes,  namely  recombination  and  translation-vibration  energy  exchange,  and  is  given 
by 

nd  nd 

Wv  =  Y^  Qt-v  s  +  X]  Wgevs  ■ 

5=1  S=1 


7.  Boundary  Conditions 

In  this  work,  we  use  slip  and  no-slip  velocity  boundary  conditions  at  the  wall.  We 
impose  either  isothermal  wall  conditions,  or  we  allow  the  internal  modes  of  the  gas  to 
’’slip”  (or  ’’jump”)  at  the  wall,  based  on  the  phenomenological  model  approach  of  Gokgen 
(1989).  The  standard  Maxwell  model  [Schaaf  &  Chambre,  (1966)  pp.  34-35]  for  velocity 
slip  and  temperature  jump  at  the  surface  of  the  body  is 


^islip 


2  -  a  .  dvt 


a  dn 


Tgljp  Tw  an  — 


2  —  ax 


wall 
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dT 


ax  (7  +  l)Pr  dn 


wall 


where  utslip  and  Tshp  are  the  tangential  velocity  and  temperature  near  the  surface  (the 
“slip”  conditions),  Twan  is  the  temperature  of  the  surface,  and  A  is  the  mean-free-path.  a 
and  ax  are  the  accommodation  coefficients  and  are  taken  to  be  0.85.  It  should  be  noted 
that  only  the  component  of  the  velocity  tangential  to  the  surface  is  allowed  to  slip  when 
velocity  slip  conditions  are  imposed.  The  component  of  the  bulk  velocity  normal  to  the 
wall  is  always  zero. 


We  also  model  the  possible  jump  of  vibrational  energy  at  the  surface.  This  is  done 
with  a  simple  extension  of  the  Maxwell  model  based  on  the  approach  of  Gokgen  (1989). 
The  vibrational  energy  jump  is 


^uslip  7;  wall 


2  a „  ^  dc„ 
av  on 


wall 


where  av  is  the  accommodation  coefficient  for  vibration  and  Xv  is  the  mean-free-path 
that  characterizes  transport  of  vibrational  energy.  Because  the  vibrational  energy  is  not 
correlated  with  the  translational  energy  modes,  we  use  the  value  of  mean-free-path  for 
the  transport  of  momentum,  namely  Xv  =  X  =  2p/pc,  where  c  =  ^8RT/tt  is  the  mean 
thermal  speed  of  the  gas.  The  value  of  the  vibrational  energy  accommodation  coefficient 
is  not  well  known,  although  Meolans  et  al.  (2002)  suggest  adiabatic  conditions  (av  =  0). 
Measurements  made  by  Black  et  al.  (1974)  give  different  values  of  av  for  different  materials. 
The  double-cone  model  used  in  this  work  is  made  of  stainless  steel,  for  which  they  measure 
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values  of  either  1.0  ±  0.2  x  10-3  or  1.2  ±  0.3  x  10-3.  In  our  calculations,  we  use  0.5  and 
1.0  x  10-3  for  av  to  assess  the  effect  of  vibrational  energy  accommodation.  Based  on  theory 
and  measurements  the  lower  value  is  more  accurate. 

In  general,  the  wall  is  assumed  to  be  non-catalytic,  which  implies  that  there  is  zero 
normal  concentration  gradient  at  the  wall  for  all  species,  and  thus  the  diffusive  flux  is  zero 
at  the  wall.  Fully  catalytic  boundary  conditions  are  also  employed  in  certain  cases,  and  we 
state  explicitly  when  these  are  used.  Assuming  that  a  wall  is  fully  catalytic  implies  that 
atomic  species  recombine  at  the  wall  to  form  molecular  oxygen  and  nitrogen. 
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The  Numerical  Method 


In  this  section  we  present  the  numerical  method  we  use  to  solve  the  coupled  partial 
differential  equations  that  describe  the  flows  of  interest.  The  conservation  laws  presented 
in  the  previous  section  were  expressed  in  their  strong  form.  However,  we  obtain  numerical 
solutions  to  this  set  of  conservation  equations  by  considering  the  weak  form  of  the  equa¬ 
tions.  We  use  a  standard  finite  volume  method  to  discretize  the  governing  equations  in 
the  approximation  space. 

We  introduce  the  finite  volume  formulation,  and  the  methods  we  use  to  evaluate  the 
inviscid  and  viscous  terms.  The  implicit  method  employed  for  time  advancement  is  also 
briefly  presented.  We  will  restrict  the  remaining  of  this  section  to  the  two-dimensional 
case  in  order  to  simplify  the  calculations,  but  it  should  be  noted  that  it  is  straightforward 
to  extend  the  numerical  method  to  three  dimensions. 

1.  The  Finite  Volume  Method 

In  this  section,  we  introduce  the  finite  volume  method  of  discretization.  We  begin 
with  the  strong  form  of  the  equations  presented  in  the  previous  section.  We  write  the 
vector  of  conserved  variables  as 

U  =  (pi , . .  • ,  pns,  pu,  pv,  Ev ,  E)t  , 

and  the  conservation  equations  can  now  take  the  compact  form 

— * 

where  F  is  the  corresponding  total  flux  vector  and  W  is  the  vector  of  source  terms 

W  =  (w1,...,wns,0,0,wv,0)T  . 

The  total  flux  can  be  split  into  the  convective  (inviscid)  and  diffusive  flux  vectors 

F  =  F!+Fv  . 

The  Cartesian  components  of  the  convective  and  diffusive  flux  vectors  are 
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/  PlUl 

P2U2 

P  —  Prisons 

~~Txy 

Qvx  “b  X/s=l  (ps^vs^s) 

\Qx  +  Qvx  -  {fu)x  +  Es=l  (P  A U«) 

P1V1 

P2V2 

Pns^ns 
~Tyx 

~ryy 

Qvy  +  Yls=l(psevsvs) 

Qy  +  Qvy  —  {ju)y  +  Xys  =  l(ps^sv5) 

In  the  finite  volume  formulation,  the  weak  form  of  the  conservation  equations  is  ob¬ 
tained  by  integrating  over  an  arbitrary  control  volume  which,  for  the  purposes  of  this 
work,  is  assumed  to  be  fixed  in  space.  Green’s  Theorem  is  applied  to  convert  the  volume 
integral  containing  the  flux  vector,  to  a  surface  integral  over  the  surface  that  encloses  12. 
The  resulting  weak  form  of  the  equation  is 

where  V  is  the  total  volume,  n  is  the  outward-pointing  normal  to  the  surface,  and  U,  W 
are  averaged  over  fh  We  discretize  physical  space  by  dividing  it  into  a  series  of  volume 
elements,  which  can  be  arbitrary  polygons.  We  do  this  by  building  a  regular  mesh  over 
the  domain  of  interest,  resulting  in  a  structured  grid.  The  grid  directions  are  defined  such 
that  i  index  increases  in  the  ^-direction  and  j  increases  in  the  positive  77-direction.  Each 
of  the  volume  elements  in  the  two-dimensional  domain  is  a  quadrilateral,  and  represents  a 
cell  in  computational  space.  Flow  variables  are  stored  at  the  cell  centers  rather  than  the 
nodes.  Quantities  at  the  cell  faces  are  computed  from  the  solution  reconstruction  based 
on  the  cell-centered  data.  The  time  rate  of  change  of  U  in  the  2-D  volume  element  i,j  is 
represented  as  a  sum  of  the  fluxes  over  the  four  cell  faces  and  the  generation  of  U  within 
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the  volume 


^(F •  h)i+ijSi+i  j  (F  ■  n)i-L ,jSi_i  j 

y  i>3  L 

+  ( F  ■  —  ( F  ■  n)ij_iSi  j_i^  +  Wij  • 

The  ±1/2  indices  on  the  flux  and  surface  area  terms  indicate  that  quantities  should  be 
evaluated  at  the  appropriate  cell  face.  This  expression  is  the  discretized  form  of  the 
conservation  equations  for  a  two-dimensional  domain.  The  extension  of  the  formulation  to 
three  dimensions,  where  the  volume  elements  are  hexahedra,  is  straight-forward  and  results 
in  an  analogous  expression.  The  axisymmetric  formulation,  however,  is  non-trivial,  and 
ultimately  yields  an  expression  similar  to  that  of  the  two-dimensional  case.  Integrating  the 
conservation  equations  over  an  axisymmetric  domain  introduces  extra  terms  that  appear 
as  sources  in  the  forcing  vector  W. 


tiki 

dt 


2.  Evaluation  of  the  Fluxes 

In  this  section  we  present  the  methods  we  use  for  evaluating  the  fluxes  that  appear 
on  the  right  hand  side  of  the  discretized  form  of  the  conservation  equations.  We  have 
decomposed  the  fluxes  into  the  convective  part  (which  consists  of  the  advection  terms 
and  the  pressure)  and  diffusive  part.  This  decomposition  is  useful  because  we  treat  the 
convective  and  diffusive  fluxes  differently.  The  inviscid  terms  make  the  system  of  equations 
hyperbolic,  while  the  viscous  terms  are  diffusive  in  nature,  and  make  the  system  parabolic 
when  combined  with  the  transient  term.  This  complicates  the  evaluation  of  the  right  hand 
side  of  the  conservation  equations.  The  elliptic  nature  of  the  diffusive  terms  allows  for  a 
simple  central  difference  scheme  to  be  used  [Hirsch  (1991)].  A  more  sophisticated  approach 
is  necessary  for  the  calculation  of  the  inviscid  terms,  in  order  to  maintain  stability  of  the 
numerical  scheme. 

In  the  computation  of  the  inviscid  terms,  we  make  use  of  the  fact  that  in  hyperbolic 
systems  of  equations,  information  propagates  along  definable  characteristic  directions  in 
space-time.  This  is  commonly  referred  to  as  upwind  biasing.  We  employ  a  modified 
form  of  the  Steger- Warming  (1981)  flux-vector  splitting,  a  method  widely  used  due  to  its 
robustness,  which  is  based  on  characteristic  theory.  We  present  the  application  of  this 
method  specifically  to  the  set  of  conservation  equations  relevant  to  this  work  (see  also 
Wright  (1997)),  but  it  should  be  noted  that  the  method  has  been  implemented  for  more 
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complex  sets  of  conservation  equations.  We  consider  the  inviscid  flux  through  a  face 


/  Pi«'  \ 
I  P2u'  I 


Fj  —  Fi  ■  n  = 
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where  sx,s'y  are  the  components  of  the  unit  vector  normal  to  the  surface  in  question,  and 
v!  is  the  component  of  the  velocity  normal  to  the  surface.  We  observe  that  the  flux  F[  is 
linear  and  homogeneous  in  the  vector  of  conserved  variables  U,  i.e., 


F'j{\U)  =  \F'j(U)  , 


where  A  is  an  arbitrary  scalar.  This  can  be  shown  if  we  express  the  inviscid  flux  Fj  in 
terms  of  the  components  of  U.  The  homogeneity  of  the  flux  allows  as  to  express  it  exactly 
in  terms  of  the  linearization 

F',(U)  =  -gJU  =  A'U  , 

where  A'  is  the  inviscid  flux  Jacobian  matrix  in  the  direction  of  h.  We  follow  the  approach 
of  Steger  and  Warming  to  split  the  inviscid  fluxes  into  positively-moving  and  negatively- 
moving  components  using  an  upwind  biasing  scheme.  The  partitioning  is  performed  by 
diagonalizing  A',  and  the  fluxes  are  split  according  to  the  sign  of  the  eigenvalues  of  the 
Jacobian.  The  diagonalization  of  is  rather  cumbersome,  but  we  can  simplify  the 
operation  by  breaking  the  Jacobian  into  components,  as 

,  _  dU  dV  dFj  8V 
~  dV  dU  dV  dU  ’ 

where  V  is  a  vector  of  primitive  variables,  introduced  for  convenience.  The  choice  of  V  is 
not  unique,  but  for  the  particular  inviscid  fluxes  it  is  convenient  to  use 

V  (pi ,  P2 ,  •  •  •  ,  pns >  j  P)  > 

where  ev  is  the  total  vibrational  energy  per  unit  mass.  The  transformation  matrices 
between  primitive  and  conservative  variables  are  denoted  by  5-1  and  S.  We  should,  in 
principle,  be  able  to  diagonalize  using  a  similarity  transformation 

dVdF'j_  i 

dU~dv~cA'  Aa'Ca'  ’ 
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where  A^'  is  a  diagonal  matrix  of  eigenvalues,  is  the  matrix  of  eigenvectors  of  A', 
and  CA>  such  that  G~^Ca>  =  I-  The  diagonalization  process  has  now  been  simplified, 
however  the  resulting  expressions  for  the  eigenvectors,  although  mathematically  correct, 
may  exhibit  singularities  and  are  not  suitable  for  use  in  a  numerical  algorithm.  Therefore, 
we  introduce  an  orthogonal  transformation  R,  which  transforms  the  momenta  from  a 
coordinate  system  locally  aligned  with  the  face  in  question,  to  the  Cartesian  system.  We 
can  now  write  the  Jacobian  matrix  as 

A'  =  S-'R-'Ca'Aa'CaR  S  , 


where  Ca  does  not  exhibit  singularities.  The  diagonal  components  of  A  a1  are  the  con¬ 
vection  speeds  of  the  characteristic  variables.  The  flux  vector  can  now  be  split  into  the 
positively  and  negatively  moving  components  defined  by 

F'+  =  S-lR-lC^K+CAR  SU  =  A'+U 
FL  =  S-'R-'Ca'A-CaR  SU  =  A'_U  , 

where  A+  and  A_  are  the  split  eigenvalue  matrices  consisting  of  the  positive  and  negative 
eigenvalues  respectively.  The  total  inviscid  flux  vector  through  the  face  in  question  can  be 
expressed  as  the  sum  of  positively  and  negatively  moving  components 

F'j  =  F’+  +  FL  . 

In  this  basic  discretization,  the  split  fluxes  at  each  cell  face  are  evaluated  using 

F+i+%,j  =  A+ijUij  , 

F~i+h>i  =  A-i+ijUi+lj  > 

where  A'±  and  U  are  evaluated  based  on  the  upwind  cell  data.  MacCormack  and  Candler 
(1989)  proposed  a  modification  to  this  method,  which  reduces  the  amount  of  numerical 
dissipation  of  the  original  method.  Their  modification  suggests  that  the  Jacobian  matrix 
should  be  evaluated  at  the  cell  face,  rather  than  using  the  upwind  cell  data.  This  approach 
works  well  in  regions  of  weak  gradients,  but  additional  dissipation  is  required  to  capture 
strong  gradients  and  shock  waves.  In  this  work,  we  use  a  hybrid  method,  in  which  we 
use  pressure  weighted  average  quantities  between  adjacent  cells  to  evaluate  the  Jacobian. 
The  method  smoothly  switches  from  modified  to  true  Steger- Warming  in  regions  of  large 
pressure  gradients.  This  flux  evaluation  method  as  originally  formulated  is  only  first  order 
accurate  in  space,  however  it  is  possible  to  achieve  higher  order  of  accuracy.  Higher  order 
accurate  schemes  can  produce  more  resolved  solutions  on  a  given  computational  mesh.  We 
use  up  to  second  order  accurate  fluxes  in  this  work.  We  obtain  higher  order  of  accuracy  by 
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extrapolating  the  conserved  variables  to  the  cell  face  based  on  the  purely  upwind  cell  data. 
We  detect  the  presence  of  shock  waves  by  identifying  pressure  jumps  greater  than  50%  of 
the  smallest  pressure  involved  in  the  flux  evaluation  stencil,  in  which  case  we  switch  the 
stencil  to  first  order  For  the  details  of  this  formulation  see  Druguet  et  al.  (2003). 


If  we  substitute  the  split  fluxes  into  the  discretized  form  of  the  conservation  equations, 
we  get  the  upwind  finite  volume  formulation  of  the  inviscid  equations 


dUij 

dt 


~  (^-  *- i  ,j  Bi~  J  ,j  Ui,j  -  Al_i+±jSi+ljUi+ Ij) 

+  (B+i,j+%Si>3+hUi’j  - 

-(B-i,j-hSh3~hUi'j  ~  B-i,j+±Si,3+hUi>3+l)]  + 


Wi 


i,3  > 


where  we  have  implicitly  defined  A'  and  B'  to  denote  the  inviscid  flux  Jacobians  in  the  £ 
and  r)  directions  respectively. 


The  elliptic  nature  of  the  viscous  fluxes  makes  their  evaluation  possible  with  a  simple 
finite  difference  scheme.  The  viscous  flux  vector  is  evaluated  at  the  required  cell  faces 
using  central  differencing.  The  viscous  fluxes  integrated  over  a  face  take  the  simple  form 


Bvi+%,jBi+\,3  ■”)*+! ,jBi+%,3 

—  (Fi  ■  ’ 


The  contribution  of  the  viscous  fluxes  is  summed  over  the  four  cell  faces  and  added  to  the 
right  hand  side  of  the  discretized  form  of  the  conservation  equations.  The  required  deriva¬ 
tives  that  appear  in  the  viscous  fluxes  are  taken  with  respect  to  the  general  curvilinear 
coordinate  system  £-77,  and  are  then  transformed  to  the  Cartesian  system.  The  calculation 
of  these  terms  can  be  found  in  Hirsch  (1991). 


3.  Time  Advancement 

In  this  section  we  discuss  how  the  solution  is  integrated  in  time  to  its  steady  state. 
We  are  interested  solely  in  the  steady  state  solution,  if  one  exists,  and  thus  time  accuracy 
is  not  important  for  our  computations.  Therefore,  we  employ  an  implicit  method,  which 
allows  for  large  time  steps  to  be  taken,  while  maintaining  stability  of  the  numerical  scheme. 
We  give  a  short  description  of  the  implicit  formulation  here.  We  first  introduce  the  explicit 
method  of  solution,  and  then  we  extend  it  to  the  implicit  formulation. 

We  discretize  the  transient  term  that  appears  in  the  conservation  equations,  in  order  to 
time-march  the  solution  toward  the  steady  state.  We  use  simple  first-order  Euler  backward 
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differencing,  in  which  the  time  derivative  becomes 


dUi 


dt 


UTtf1-U*.  A  U?1- 


and  the  superscript  denotes  the  time  level  in  the  temporal  discretization.  This  expression 
is  substituted  directly  into  the  discretized  form  of  the  conservation  equations,  and  the 
solution  is  marched  forward  in  time  by  discrete  time  steps.  The  explicit  formulation  of  the 
inviscid  problem  is 


At  r 

■v~ 


(A+ i+  i  ,j  Si+  3  >J  Ui>j  A'+i-\  ,j  d  Ui_1  ’j  ) 

(A'-i-ijSi-ijUij  -  A'_.+  i.Si+ijUi+ij) 


+  (B+i,j+%Sh3+hUi’3  1) 

+  AtWTj  , 


where  A 17” ■  is  the  explicit  change  in  the  solution  vector  (residual)  at  cell  i,j  at  time  step 
n.  The  solution  is  updated  to  the  next  time  level  using 

v r/1  =  v?j + w?j . 


We  call  this  an  explicit  method  because  the  solution  in  the  next  time  level  n+ 1,  is  computed 
explicitly  by  evaluating  the  right  hand  side  at  the  current  time  level  n.  The  solution  at 
any  point  at  time  level  n  +  1  is  independent  of  the  solution  at  other  points  at  the  same 
time  level. 

With  a  fully  implicit  method,  the  time-step  limitations  inherent  to  explicit  solvers 
can  be  overcome.  The  use  of  an  implicit  solver  can  allow,  in  general,  for  time  steps  much 
larger  than  the  stable  explicit  time  step.  In  an  implicit  method,  the  solution  at  any  point 
depends  on  the  solution  of  all  other  points  at  the  time  level  n  +  1.  Therefore,  in  order  to 
employ  implicit  time  advancement,  it  is  necessary  to  evaluate  the  right  hand  side  of  the 
conservation  equations  at  time  level  n  +  1.  We  express  the  fluxes  that  appear  on  the  right 
hand  side  in  terms  of  the  linearization 


jpf  71+1 


dF[ 


+  dU 
+  A'n8U 


{Un+ 1  -  Un)  +  0(At2) 


where  we  have  defined  the  vector  6Un  —  Un+1  —  Un  and  we  have  frozen  the  variation  of 
the  inviscid  flux  Jacobian  A'n  between  time  steps.  We  apply  this  linearization  to  the  split 
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fluxes  and  source  vector  to  obtain  the  fully  implicit  upwind  formulation  of  the  inviscid 
problem 

5U"i  +  Sr  ^  jWi-u) 

-  ALiHtjS,+hs5UWJ) 
+K«+4sy+Jw«  -  Kij-iSij-iWij-i) 

-  S'-y+jSy+}Wu+,)]” 

-At  CTJU^i  =  At/£  , 

where  C!jb  is  the  Jacobian  of  W  with  respect  to  U.  In  this  implicit  equation,  A £/£.  is  the 
explicit  residual  as  defined  previously. 

We  include  the  implicit  approximation  of  the  viscous  terms  in  the  implicit  equation 
by  applying  the  linearization  to  the  viscous  fluxes.  We  define  5F„n  and  SG'vn  such  that 


p/  n+1 


F'+SF' 


G'vn+1  =  G'vn  +  6G'vn 


Therefore,  we  only  need  to  find  expressions  for  5F^  and  5G'V.  We  note  that  and  G'v 
are  functions  of  U  and  its  spatial  derivatives.  The  derivatives  are  taken  with  respect  to 
the  general  curvilinear  coordinate  system,  and  then  transformed  to  the  Cartesian  system 
using  the  operators 

d  _  d£  d  dr)  d  d  _  d£  d  dr)  d 

dx  dx  3£  dx  dr)  ’  dy  dy  <9£  +  dy  dr) 

We  apply  the  thin-layer  assumption  to  the  derivatives  in  the  implicit  viscous  terms.  In 
this  way,  we  selectively  retain  derivatives  depending  on  which  viscous  flux  we  wish  to 
approximate.  For  example,  in  the  approximation  of  5G'V  we  compute  the  derivatives  by 

d  ^  drj  d  d  dr)  d 

dx  dx  dr)  ’  dy  ~  dy  dr) 

Similarly,  derivatives  with  respect  to  £  are  retained  when  approximating  8F'V.  If  we  assume 
that  the  terms  not  involving  derivatives  are  locally  constant,  we  can  write  F'v  and  G'v  as  a 
function  of  the  derivatives  of  the  flow  variables  with  respect  to  £  or  77  only,  as 

5F'V  =  M^(6V)  ,  SG'V  =  Mv-^-(5V)  , 

and  we  introduce  for  convenience  the  vector  of  non-conserved  variables 


V  =  (ci,c2, . . .  ,cns,u,v, e„,T)T  . 
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The  matrices  M j  and  M, 7  can  be  found  in  Candler  (1988).  We  ultimately  express  the 
implicit  viscous  fluxes  in  terms  of  5U  by  mapping  from  5V  to  6U  using  the  Jacobian 
N  =  ^ ,  and  we  obtain  the  final  form 


5F'V  =  M^{N5U)  ,  5G'V  =  MV-^(N5U)  , 

We  should  emphasize  here  that  all  of  the  viscous  derivatives  are  retained  in  the  computation 
of  F^n  and  G'vn ,  and  thus  all  of  the  viscous  terms  influence  the  converged  solution.  With 
this  treatment  of  the  viscous  terms,  the  implicit  equation  for  the  inviscid  problem  away 
from  solid  boundaries  will  be  unchanged  if  we  simply  replace  the  inviscid  Jacobians  A'  and 
B'  with  A  and  B,  where. 

A+  =  A+  -  M^N  A-  =  A-  +  M^N 
B+  =  B+-  MVN  +  Mr,N 

The  implicit  formulation  involves  the  solution  of  a  linear  system  of  equations,  which 
has  a  strong  diagonal  dominance  due  to  the  nature  of  the  Steger- Warming  Jacobians.  In 
principle,  we  can  solve  exactly  for  the  by  a  direct  inversion  of  the  operator. 

4.  The  Data-Parallel  Line-Relaxation 

In  this  section  we  briefly  present  the  method  of  solution  of  the  implicit  equation.  We 
use  the  Data-Parallel  Line  Relaxation  method  of  Wright  et  al.  (1998).  This  method  was 
developed  to  take  advantage  of  parallel  computing  machines,  and  appears  to  be  superior 
to  other  implicit  methods  [Wright  (1997)].  Only  a  brief  description  of  this  method  is  given 
here.  Details  on  the  performance  of  the  method  can  be  found  in  Wright  (1997)  and  Wright 
et  al.  (1998). 

The  standard  implicit  upwind  formulation  of  the  Navier-Stokes  equations  for  a  reacting 
gas  is  essentially  a  linear  system  of  equations.  The  solution  of  such  system  can,  in  principle, 
be  obtained  directly  by  inversion  of  a  large  block  banded  matrix.  This  becomes  very 
computationally  intensive,  as  such  an  inversion  must  be  performed  at  every  time  step  until 
convergence  is  achieved.  Therefore,  most  implicit  methods  seek  to  make  simplifications  in 
order  to  make  the  solution  of  the  problem  less  expensive.  Such  simplifications  are  possible 
only  on  the  left-hand  side  of  the  equations,  which  resembles  the  numerics.  The  right-hand 
side  of  the  equation  (residual)  should  be  unaltered. 

It  is  convenient  to  introduce  some  additional  notation  in  order  to  group  terms  on  the 
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left-hand  side  of  the  implicit  equation.  We  let 


Ai,j  — 


A  tI+  iA  +^+5  jSi+i>3 

-A,n_.  ,  .Si_ i, 

*-5>J  *  2’J 

-  VijCT. 

0/n  Q 

_  ~n,n  c 

Ui<j  ~  **  5 

Di,  =  A'li+hjSi+i>j 


Ei.<=A'l;  , 


hi 


We  can  simplify  the  operator  if  we  assume  that  the  physical  problem  is  strongly 
coupled  in  the  rj  direction.  Then,  we  can  move  some  of  the  off-diagonal  terms  to  the 
right-hand  side  of  the  implicit  equation,  keeping  only  constant  i  terms  on  the  left-hand 
side.  In  this  way,  it  is  possible  to  solve  for  all  j  points  at  each  i  location  as  a  series  of  fully 
coupled  block  tri-diagonal  systems  aligned  in  the  rj  direction.  We  relax  the  off-diagonal 
terms  by  updating  the  unknowns  on  the  right-hand  side  of  the  equation  with  the  SUij  we 
have  just  computed.  We  do  not  need  to  obtain  an  exact  solution  to  the  implicit  equation, 
but  rather  to  maintain  stability,  and  thus  only  a  few  iterations  axe  required  per  time  step. 
The  solution  procedure  is  as  follows.  First,  the  implicit  terms  on  the  right-hand  side 
are  neglected  and  the  resulting  block  tri-diagonal  system  is  factored  and  solved  for  SU^ 
according  to: 

Bi,M%  +  AMy  -  Cijtuft =  A U?J 

where  the  superscript  denotes  the  sub-iteration  number  (zero  denotes  the  initialization). 
Then,  a  series  of  kmax  relaxation  steps  is  performed,  where  we  solve 

+  hMj  -  +  AC'S  . 

and  the  superscript  k  —  1  indicates  that  data  from  the  previous  relaxation  step  is  used. 
Finally,  the  solution  is  obtained  by 


suff1  -  suQmax) . 


There  are  two  advantages  to  this  method  of  solution.  First,  the  method  can  be  imple¬ 
mented  efficiently  on  a  distributed  memory  parallel  computer.  This  is  because  we  can  split 
the  problem  among  processors  in  the  i  direction,  effectively  assigning  a  series  of  i-columns 
to  each  processor.  Therefore,  each  relaxation  step  can  be  performed  simultaneously  as 
there  are  no  data  dependencies,  while  the  boundary  data  can  be  efficiently  transfered. 
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Secondly,  there  is  no  bias  in  the  solution  procedure,  which  has  been  shown  to  impose  a 
bias  in  the  solution  itself. 

5.  Boundary  Conditions 

In  this  section  we  describe  the  boundary  conditions  used  for  both  the  explicit  residual 
and  the  implicit  operator,  and  how  they  are  implemented  in  the  numerical  algorithm. 

In  the  finite  volume  method,  explicit  boundary  conditions  are  imposed  using  “dummy” 
cells.  An  extra  layer  of  cells  is  included  around  the  grid  that  represents  physical  space,  in 
which  variables  are  stored  in  order  to  impose  the  appropriate  boundary  conditions  for  the 
problem.  The  dummy  cells  are  constructed  in  such  a  way  that  derivatives  taken  in  terms 
of  the  curvilinear  body-fitted  coordinate  system  are  consistent  with  the  type  of  boundary 
condition  of  the  viscous  fluxes.  Because  we  split  the  fluxes  into  the  viscous  and  inviscid 
parts,  we  treat  the  boundary  conditions  separately. 

In  the  explicit  formulation,  we  simulate  no-slip  conditions  by  setting  the  velocity  at 
the  wall  to  zero,  and  extrapolating  the  velocity  components  at  the  boundary  cell  i,  1  based 
on  the  values  at  the  cell  i, 2.  We  simulate  isothermal  wall  conditions  in  a  similar  fashion. 
The  slip  conditions  require  that  we  have  a  model  for  how  the  velocity  component  tangent 
to  the  wall  is  computed  based  on  the  interior  values.  The  details  of  the  model  employed 
were  presented  in  the  previous  section.  Once  the  slip  velocity  has  been  determined,  we 
extrapolate  values  to  the  dummy  cells,  as  we  did  for  the  no-slip  case.  We  impose  internal 
energy  jumps  at  the  wall  in  exactly  the  same  fashion.  The  non-catalicity  of  the  wall  is 
imposed  by  setting  the  values  of  the  species  densities  at  the  dummy  cells  equal  to  the 
interior  values,  effectively  imposing  a  zero  mass-fraction  gradient.  Similarly,  fully  catalytic 
boundary  conditions  are  imposed  be  extrapolating  mass-fractions  to  the  dummy  cell  as¬ 
suming  the  mass-fractions  of  atomic  species  is  zero  at  the  wall.  Subsequently,  the  mass 
flux  of  polyatomic  species  at  the  wall  is  corrected  such  that  we  have  elemental  conservation 
during  this  process. 

The  implicit  treatment  of  the  boundary  conditions  is  analogous  to  that  of  the  explicit 
method.  In  the  implicit  formulation,  however,  it  is  not  the  variables  at  dummy  cells  that 
are  set  to  a  certain  value,  but  their  effect  on  the  implicit  operator  which  we  construct. 
Specifically,  the  boundary  conditions  are  folded  into  the  implicit  operator  by  making  a 
modification  to  the  Aij  matrix  where  appropriate.  For  example,  the  inviscid  boundary 
conditions  at  the  wall  are  included  by  adding  the  matrix  Ci^Ew  to  A,^-  where  Cij  was 
defined  previously,  and  Ew  is  an  operator  that  reflects  the  momentum  normal  to  the  wall. 
Similarly,  the  approximate  implicit  viscous  Jacobian  Mn  is  modified  at  solid  boundaries. 
The  details  of  the  implicit  treatments  of  boundaries  are  described  in  Candler  (1988)  and 
G  ok  gen  (1989).  Note  that  the  viscous  wall  implicit  boundary  conditions  are  folded  into 
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the  block  matrix  Ai$  and  should  not  be  included  in  Cit 2  • 

6.  Other  Flux  Evaluation  Methods 

The  inviscid  flux  evaluation  method  discussed  earlier  is  one  of  many  methods  suitable 
for  gasdynamic  simulations.  There  are  several  methods  used  widely  in  numerical  algo¬ 
rithms,  and  these  are  of  variable  degree  of  sophistication  and  complexity.  Different  flux 
evaluation  methods  exhibit  different  properties  in  terms  of  accuracy  and  computational 
cost.  Most  importantly,  the  amount  of  numerical  dissipation  that  is  associated  with  a  nu¬ 
merical  scheme  depends  largely  on  the  choice  of  flux  evaluation  and  solution  reconstruction 
methods.  In  general,  a  scheme  with  little  numerical  dissipation  produces  more  accurate 
results  on  a  given  computational  mesh. 

It  is  beyond  the  scope  of  this  work  to  make  comparisons  between  different  flux  eval¬ 
uations  methods.  However,  it  is  important  to  assess  how  well  the  numerical  method  of 
choice  performs  in  comparison  with  other  methods.  Therefore,  a  later  section  discusses 
in  more  detail  results  of  simulations  made  using  several  other  methods  of  evaluating  the 
inviscid  fluxes.  These  flux  evaluation  methods  are  not  covered  in  great  detail,  and  are  only 
presented  briefly. 

7.  Second  Order  Spatial  Accuracy 

Higher  order  accurate  inviscid  fluxes  can  be  obtained  through  appropriate  reconstruc¬ 
tion  of  the  solution  using  cell-centered  data.  Based  on  the  choice  of  the  reconstruction 
method,  the  accuracy  of  the  inviscid  fluxes  may  differ.  In  this  section,  we  briefly  present 
the  methods  used  in  this  work. 

The  flux  evaluation  method  discussed  previously  is  first  order  accurate  when  data 
from  adjacent  cells  are  used  to  calculate  the  flux.  The  method  can  become  second  order 
accurate  with  a  simple  modification.  Consider  the  Steger- Warming  flux  for  the  inviscid 
equations  in  one  spatial  dimension: 

F  =  F+(Ui)  +  F-(Ui+1)  . 

In  this  form,  the  split  fluxes  use  upwind  data  from  adjacent  cells.  A  higher  order  numerical 
flux  H  may  be  obtained  if  we  replace  the  arguments  of  F+  and  F_  with  a  left  state  UL 
and  a  right  state  UR  as: 

H  =  F+(Ul)  +  F-(Ur) 

=  A+Ul+A-Ur  , 

where  UL  and  UR  are  extrapolated  conserved  quantities  based  on  pure  upwind  data.  It 
should  be  noted  that  this  simple  extrapolation  is  not  suitable  for  regions  of  the  flow  that 
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exhibits  large  gradients  in  f7,  such  as  shock  waves.  Therefore,  we  need  to  detect  the 
presence  of  a  shock,  and  use  no  extrapolation  in  that  case.  Our  implementation  uses 
pressure  jumps  across  cell  faces  in  order  to  detect  the  presence  of  shocks,  as  discussed 
earlier. 

The  MUSCL  variable  extrapolation  method  may  also  be  used  to  construct  higher 
order  Steger- Warming  and  modified  Steger- Warming  fluxes.  In  this  case,  the  left  and  right 
states  UL,  UR  are  constructed  based  on  non-linear  functions  of  the  ratio  of  adjacent  slopes 
of  U  computed  based  on  cell-centered  data.  Details  on  slope  limiters  and  their  use  in  the 
MUSCL  variable  extrapolation  can  be  found  in  Hirsch  (1991)  and  Laney  (1998). 
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Low  Enthalpy  Double-Cone  Experiments 

In  this  section  we  present  simulations  of  double-cone  experiments  performed  at  low 
enthalpy  conditions,  and  we  show  blind  comparisons  between  the  experimental  measure¬ 
ments  and  numerical  predictions.  The  term  blind  implies  that  there  was  no  prior  knowledge 
of  the  measurements  obtained  experimentally  when  the  simulations  were  performed.  The 
experimentalists  had  not  released  the  results  before  the  simulations  were  completed,  and 
only  the  free-stream  conditions  and  the  model  geometry  specification  were  provided.  In 
this  section,  we  also  briefly  describe  how  the  experiments  were  performed,  and  give  a  short 
description  of  the  experimental  facility. 

The  experiments  were  part  of  a  blind  validation  exercise,  that  was  conducted  under 
the  auspices  of  the  NATO  Research  and  Technology  Organization  (RTO)  Working  Group 
10  [Holden  et  al.  (2001)  &  Harvey  et  al.  (2001)].  This  exercise  was  part  of  a  large 
coordinated  effort  by  the  working  group  to  assess  the  ability  of  popular  CFD  codes  to 
reproduce  hypersonic  flows  that  closely  resemble  what  is  encountered  during  actual  flight. 
This  effort  is  based  on  a  building  block  approach,  where  different  aspects  of  the  simulation 
are  tested  and  validated  progressively.  In  particular,  the  numerical  methods  and  the  ability 
of  simulations  to  accurately  reproduce  complex  shock  interaction  phenomena  in  essentially 
perfect-gas  flows  were  investigated  first.  Then,  depending  on  the  outcome  of  the  first 
step,  more  complex  hypersonic  flow  phenomena  requiring  modeling  of  physical  processes 
were  examined.  One  of  the  objectives  of  the  working  group  was  to  build  a  database  of 
experimental  measurements  for  flows  that  exhibit  various  physical  phenomena  -  such  as 
high  enthalpy  effects  and  turbulence  -  along  with  benchmark  numerical  predictions  by 
continuum  and  particle-based  methods.  The  low  enthalpy  experiments  presented  here 
are  part  of  validation  efforts  focusing  on  the  numerical  aspect  of  simulations,  and  aim  to 
address  issues  of  numerical  accuracy  versus  computational  cost. 

1.  Low  Enthalpy  LENS  Experiments 

The  experiments  presented  in  this  section  are  of  hypersonic  flows  of  nitrogen  over 
a  25°-55°  double-cone  at  free-stream  specific  enthalpies  of  about  3.5  MJ/kg  and  nominal 
Mach  numbers  of  about  10.  The  experiments  were  performed  in  the  Large  Energy  National 
Shock  tunnel  (LENS)  facility  at  Calspan  -  University  at  Buffalo  Research  Center  (CUBRC), 
using  the  LENS  Leg  I  reflected  shock  tunnel.  The  tunnel  consists  of  a  shock  tube  connected 
to  a  hypersonic  nozzle  leading  to  the  test-section.  A  diaphragm  divides  the  shock  tube  into 
two  compartments.  The  first  compartment  contains  the  driver  gas,  which  can  be  helium 
or  hydrogen.  The  second  compartment  contains  the  test  gas.  The  test  gas  is  shock  heated 
twice;  first  by  the  normal  shock  that  is  formed  when  the  diaphragm  raptures,  and  then  as 
the  shock  reflects  from  the  end  of  the  tube.  The  stagnant  gas  behind  the  reflected  shock 
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breaks  the  secondary  diaphragm,  escapes  through  the  nozzle  throat  and  expands  in  the 
nozzle  to  hypersonic  conditions.  Measurements  are  made  throughout  the  experiment,  both 
inside  the  shock  tube  and  at  the  exit-plane  of  the  nozzle.  Based  on  these  measurements, 
the  free-stream  conditions  are  computed  by  expanding  the  measured  reservoir  (stagnation) 
conditions  assuming  quasi  one-dimensional  flow. 

These  experiments  were  designed  to  result  in  low  free-stream  Reynolds  numbers,  in 
order  to  exhibit  fully  laminar  flow  over  the  double-cone  model.  Additionally,  the  free- 
stream  enthalpy  is  low  enough  such  that  there  is  no  dissociation.  We  present  simulations 
for  two  experimental  conditions  of  several  available  -  Run  28  and  Run  35.  These  two 
sets  of  conditions  are  of  nitrogen  flows  at  similar  Reynolds  numbers,  but  they  differ  in 
the  free-stream  Mach  numbers.  The  free-stream  and  wall  conditions  for  these  cases  are 
tabulated  in  Table  1. 


Case 

Run  28 

Run  35 

Moo 

9.59 

11.30 

Rem  (103/m) 

130.9 

133.3 

Too  (K) 

185.56 

138.89 

Twall  (K) 

293.33 

296.11 

Poo  (10-3  kg/m3) 

0.6545 

0.5515 

Uoo  (m/s) 

2663.9 

2712.7 

Table.  1  Free-stream  and  wall  conditions  for  the  25°-55°  sharp  double-cone  LENS 
experiments. 

The  available  data  for  making  comparisons  between  numerical  results  and  the  ex¬ 
periments  are  surface  quantities  and  Schlieren  images.  The  surface  quantities  are  surface 
pressure  and  heat  transfer  rate,  measured  at  several  locations  along  the  model.  There  are 
approximately  25  pressure  measurements  and  40  heat  transfer  rate  measurements  for  each 
of  the  two  cases.  This  large  number  of  measurements  makes  these  experiments  the  high¬ 
est  resolution  data  obtained  over  double-cone  flows.  It  should  be  noted  that  such  a  large 
number  of  measurements  is  desirable  because  it  enables  us  to  get  an  accurate  estimate  of 
the  separation  zone  size. 

2.  Computations  of  the  Experiments 

In  this  section  we  present  results  from  numerical  simulations  of  the  LENS  low  enthalpy 
experiments  and  we  make  direct  comparisons  with  the  experimental  data.  We  emphasize 
that  the  results  presented  in  this  section  are  compared  against  the  experimental  measure¬ 
ments  in  a  blind  manner,  as  there  had  been  no  prior  knowledge  of  the  experimental  results 
when  the  simulations  were  performed.  The  experimental  data  were  released  after  comple¬ 
tion  of  the  simulations,  and  after  the  first  comparisons  were  made  by  a  third  party  (see 
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Figure  2.  Computed  Mach  number  contours  of  the  flow  generated  by  a  25°-55° 
double-cone  geometry  at  Run  28  conditions. 

Harvey  et  al.  (2001)). 

We  simulated  the  flow  assuming  that  the  ffee-stream  was  in  equilibrium,  and  we 
allowed  for  vibrational  excitation  and  nonequilibrium  nitrogen  dissociation  based  on  the 
Park  model.  We  observed  a  very  small  degree  of  vibrational  excitation,  which  accounted 
for  a  small  fraction  of  the  energy.  There  was  essentially  no  dissociation  in  the  flow.  We  also 
performed  a  grid  convergence  study  and  found  that  we  obtained  grid-independent  results 
with  a  mesh  of  at  least  1024  x  512  points  in  the  streamwise  and  wall-normal  directions. 
Therefore,  we  used  a  1024  x  512  point  mesh  for  these  simulations.  The  simulations  were 
performed  on  a  cluster  of  DEC/Compaq  high-performance  XP-1000  Alpha  workstations 
connected  with  a  Myrinet  interconnect.  We  obtained  perfect  scaling  on  the  finest  mesh  up 
to  the  16  available  processors.  These  simulations  required  about  24  hours  of  simulation 
time. 

Fig.  2  shows  computed  Mach  number  contours  under  Run  28  conditions.  The  flow  in 
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Figure  4.  Computed  surface  pressure  and  experimental  measurements  plotted 
versus  axial  distance  on  the  25°-55°  double-cone  geometry  at  Run  28  conditions. 

this  figure  is  circularly  symmetric  about  a  horizontal  axis  passing  through  the  tip  of  the 
geometry.  The  darkest  contour  shade  indicates  regions  of  the  domain  where  the  flow  is 
subsonic.  The  oblique  shock  originating  from  the  sharp  tip  is  visible.  The  shock  generated 
by  the  separated  region  of  the  flow  interacts  with  the  bow  shock,  and  the  triple  point  is 
clearly  defined.  Also  shown  is  the  subsonic  region  behind  the  strong  bow  shock,  as  well 
as  the  supersonic  jet.  A  small  kink  is  visible  along  the  bow  shock.  The  peak  temperature 
is  observed  behind  the  bow  shock  near  the  triple  point,  where  the  flow  goes  through, 
essentially,  a  normal  shock.  The  free  shear-layer  that  forms  at  the  edge  of  the  separation 
is  highlighted  by  the  change  in  the  Mach  number  across  the  flow  direction. 

Fig.  3  shows  an  experimental  and  numerical  Schlieren  image  of  the  flow  over  the  double 
cone  for  Run  28.  Because  of  the  large  size  of  the  model,  the  experimental  Schlieren  image 
shows  only  the  upper  part  of  the  double  cone.  The  numerical  Schlieren  was  constructed  by 
simply  visualizing  the  second  gradient  of  the  density  field  along  the  axisymmetric  domain. 
The  shock  structure  is  visible  in  the  experimental  image.  Some  detail  of  the  shock  structure 
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Figure  5.  Computed  heat  transfer  rate  and  experimental  measurements  plotted 
versus  axial  distance  on  the  25°-55°  double-cone  geometry  at  Run  28  conditions. 

at  the  point  of  interaction  is  also  visible.  Although  the  numerical  image  looks  qualitatively 
similar  to  the  experimental  image,  it  appears  that  the  triple  point  is  shifted  downstream 
in  the  computed  result.  Furthermore,  the  shock  shape  is  slightly  different  in  the  vicinity 
of  the  triple  point. 

Direct  comparisons  between  the  numerical  predictions  and  the  measurements  are  made 
by  plotting  surface  quantities.  Fig.  4  shows  computed  surface  pressure  and  experimental 
measurements  plotted  versus  the  axial  distance  on  the  model.  We  see  a  region  of  constant 
pressure  along  the  first  cone  where  the  flow  is  attached.  The  sharp  increase  in  surface 
pressure  indicates  separation.  In  this  case,  the  flow  separates  at  about  half  the  length  of 
the  first  cone.  There  is  a  pressure  plateau  inside  the  separation  zone,  which  is  followed  by  a 
very  large  increase  at  the  point  where  the  transmitted  shock  impinges  on  the  surface  of  the 
second  cone.  The  fluctuations  in  surface  pressure  downstream  of  the  point  of  interaction 
are  due  to  the  compression  and  expansion  of  the  jet  near  the  surface.  The  constant  pressure 
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Figure  6.  Computed  Mach  number  contours  of  the  flow  generated  by  a  25°-55° 
double-cone  geometry  at  Run  35  conditions. 

region  along  the  first  cone  is  predicted  by  the  simulation  to  within  5%  of  the  measurements, 
This  is  well  within  the  experimental  uncertainty.  The  pressure  is  also  matched  inside  the 
separation  zone,  and  the  peak  pressure  at  the  point  of  interaction  agrees  well  with  the 
experimental  data.  The  general  trend  and  level  of  pressure  fluctuation  after  re-attachment 
are  also  reproduced,  but  are  not  exactly  matched.  This  is  because  the  location  of  the  point 
of  interaction  is  not  predicted  exactly.  We  observe  similar  results  with  the  heat  transfer 
rate  data,  which  is  shown  in  Fig.  5.  The  simulation  reproduces  the  level  of  heat  transfer 
rate  everywhere,  with  the  exception  of  the  attached  region  of  the  flow  near  the  tip.  In  this 
region,  the  heat  transfer  rate  is  over-predicted  by  the  simulation  by  about  20%.  Also,  it  is 
not  clear  whether  the  separation  point  is  predicted  exactly,  however  the  computed  result 
is  at  most  within  4  mm  of  that  observed  based  on  the  pressure  and  heat  transfer  rate  data. 

The  results  of  the  simulations  under  Run  35  conditions  are  similar  to  the  results  of 
Run  28.  Computed  Mach  number  contours  for  Run  35  are  shown  in  Fig.  6.  The  basic  flow  is 
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Figure  8.  Computed  surface  pressure  and  experimental  measurements  plotted 
versus  axial  length  of  the  25°-55°  double-cone  geometry  at  Run  35  conditions. 

similar  to  Run  28,  however  the  separated  region  is  significantly  smaller  for  this  case.  The 
darkest  contour  level  indicates  subsonic  flow,  and  this  is  seen  behind  the  strong  bow  shock 
and  in  the  separation  zone.  The  shape  of  the  supersonic  jet  and  the  kink  in  the  bow  shock 
are  also  visible.  Fig.  7  shows  the  corresponding  experimental  and  numerical  Schlieren 
images  of  the  flow  over  the  double  cone  for  Run  35.  The  shock  shapes  in  these  images 
look  remarkably  similar.  The  location  of  the  kink  in  the  bow  shock  is  predicted  very  well. 
Also,  the  shock  structure  near  the  triple  point  and  the  transmitted  shock  look  the  same, 
to  the  level  of  detail  and  the  extent  of  comparison  that  can  be  made  using  this  image. 

Fig.  8  plots  computed  surface  pressure  for  Run  35  together  with  experimental  mea¬ 
surements  versus  axial  distance.  The  pressure  in  the  attached  region  along  the  first  cone 
is  slightly  under-predicted  by  the  simulation,  but  is  within  the  5%  quoted  experimental 
uncertainty.  However,  the  level  of  pressure  inside  the  separation  zone  is  not  reproduced. 
The  difference  in  this  region  is  approximately  10%  of  the  measured  values.  It  is  not  clear 
whether  the  peak  pressure  is  over-predicted,  or  the  large  discrepancy  observed  is  because 
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Figure  9.  Computed  heat  transfer  rate  and  experimental  measurements  plotted 
versus  axial  length  of  the  25°-55°  double-cone  geometry  at  Run  35  conditions. 

of  not  having  predicted  correctly  the  location  of  the  interaction  point.  The  pressure  dis¬ 
tribution  is  predicted  reasonably  well  along  the  second  cone  after  reattachment.  Also,  the 
separation  zone  size  is  over-predicted  by  about  2  mm. 

The  heat  transfer  rate  for  Run  35  is  plotted  in  Fig.  9.  We  observe  good  agreement 
with  the  measurements  made  inside  the  separation  zone  and  after  reattachment,  but  the 
heat  transfer  rate  is  substantially  over-predicted  along  the  first  cone.  This  difference  is 
about  20%,  and  this  is  consistent  with  Run  28.  Also,  the  location  of  the  peak  in  heat 
transfer  rate  is  not  predicted  exactly. 

3.  Remarks  on  Blind  Comparisons 

The  simulations  of  Run  28  and  Run  35  from  the  low  enthalpy  LENS  experiments 
presented  in  the  previous  section  showed  good  overall  agreement  with  the  experimental 
measurements.  The  numerical  predictions  were  obtained  using  a  nonequilibrium  formula- 
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tion,  which  allowed  for  vibrational  energy  excitation  and  nitrogen  dissociation.  However, 
because  there  was  essentially  no  dissociation  and  there  was  little  vibrational  energy  exci¬ 
tation,  the  present  results  can  be  reproduced  using  a  perfect-gas  simulation. 

The  simulations  required  a  relatively  large  number  of  grid  points  in  order  to  achieve 
grid-convergence.  A  carefully  constructed  grid  of  at  least  1024  x  512  points  produced 
grid-independent  results.  This  is  consistent  with  earlier  work  on  flows  over  double-cones 
by  Wright  et  al.  (1999).  The  study  by  Gaitonde  et  al.  (2002)  shows  similar  grid  conver¬ 
gence.  They  also  conducted  a  careful  time  convergence  study,  and  showed  that  the  solution 
does  not  converge  until  at  least  100  characteristic  flow  times  have  been  established.  The 
simulations  presented  here  are  consistent  with  their  findings,  and  we  simulate  our  flows  to 
a  minimum  of  150  flow  times.  It  should  be  noted  that  such  simulations  on  a  grid  of  the 
size  required  for  grid-convergence  are  very  expensive. 

The  simulations  assume  that  the  flow  is  laminar,  which  is  consistent  with  the  exper¬ 
iments.  If  the  flow  were  transitional  or  turbulent,  the  simulations  would  under-predict 
the  heat  transfer  rate  measurements.  Also,  we  expect  the  flow  to  be  laminar  because 
the  Reynolds  number  based  on  the  free-stream  conditions  and  model  diameter  is  at  most 
40,900  for  the  cases  considered.  Additionally,  the  the  free  shear-layer  that  forms  at  the 
edge  of  the  separation  zone  -  and  is  the  most  unstable  region  in  the  flow  -  remains  steady. 
The  shear-layer  Reynolds  number  computed  based  on  the  Birch  and  Keyes  (1972)  defini¬ 
tion  for  free  shear-layers  formed  by  shock  interactions,  is  under  20,000,  and  is  expected  to 
be  stable. 

The  blind  comparisons  made  in  this  work  illustrate  that  numerical  methods  are  capa¬ 
ble  of  reproducing  double-cone  flows,  and  this  is  done  within  a  reasonable  amount  of  time. 
However,  it  is  necessary  to  use  an  efficient  implicit  solver.  The  surface  pressure  predictions 
in  regions  where  the  pressure  is  constant  are  at  most  within  10%  of  the  measurements,  and 
the  heat  transfer  rate  predictions  exhibit  a  similar  trend.  However,  the  heat  transfer  rate 
along  the  first  cone  is  overpredicted  in  both  cases  by  about  20%.  These  results  are  con¬ 
sistent  with  the  findings  of  other  researchers  who  have  simulated  these  flows  [Gaitonde  et 
al.  (2002),  Gnoffo  (2001),  &  Roy  et  al.  (2001)].  In  the  following  sections  we  perform  more 
detailed  simulations  of  the  entire  experiment,  to  determine  the  source  of  the  discrepancies. 
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Effects  of  Numerics  on  Double-Cone  Flows 

In  this  section,  we  address  the  numerical  aspects  of  the  simulation  of  double-cone  flows. 
The  main  issue  that  is  addressed  concerns  the  effect  of  the  numerics  on  the  interaction 
between  the  shocks  and  the  separation  zone.  We  show  that  the  shock  interaction  and  the 
separation  zone  are  very  sensitive  to  the  numerical  dissipation.  The  numerical  dissipation 
depends  on  the  accuracy  of  the  numerical  scheme  used  in  the  discretization  of  the  equations, 
and  on  the  refinement  of  the  mesh.  For  the  purposes  of  this  analysis  we  examine  only  one 
case,  Run  35,  from  the  LENS  experiments  presented  in  the  previous  section.  This  work 
was  done  in  collaboration  with  Dr.  Marie-Claude  Druguet  of  IUSTI  -  Ecole  Polytechnique 
Universitaire  de  Marseille. 

As  discussed  in  the  first  section,  the  influence  of  the  shock  angles  on  the  recirculation 
zone  is  very  important.  In  the  context  of  numerical  simulation,  this  translates  to  being 
able  to  predict  correctly  the  shock  locations  and  the  location  of  the  triple  point.  Thus,  it 
is  critical  to  resolve  the  initial  boundary  layer  growth  near  the  model  leading  edge,  and 
also  to  resolve  the  finer  flow  structure  at  the  point  of  interaction.  Therefore,  several  levels 
of  grid  refinement  are  used  until  grid  convergence  is  achieved.  We  show  that  grids  of  1/2 
million  points  give  accurate  results. 

To  investigate  the  effect  of  the  numerical  dissipation  of  the  schemes,  several  flux 
evaluation  methods  are  used.  We  compare  the  results  obtained  with  different  schemes 
by  displaying  physical  characteristics  of  the  flow  and  by  monitoring  the  evolution  of  the 
surface  quantities.  We  show  that  the  computed  results  are  very  sensitive  to  the  numerical 
flux  evaluation  method  used.  We  find  that  when  the  grid  is  fine  enough,  all  the  flux 
evaluation  schemes  give  the  same  results.  However,  this  may  require  a  very  large  number 
of  points  for  the  most  dissipative  schemes.  The  less  dissipative  schemes  give  the  final,  grid 
converged,  result  on  a  coarser  mesh,  and  therefore  are  less  costly.  We  also  show  that  since 
a  large  number  of  grid  points  is  needed,  it  is  necessary  to  employ  an  implicit  solver  for 
these  computations. 

1.  Grid  Size  Requirements 

We  performed  a  grid  convergence  study  on  the  double-cone  flow  under  Run  35  condi¬ 
tions.  At  each  level  of  grid  refinement,  the  number  of  points  in  each  directions  was  doubled. 
For  these  simulations,  we  used  the  second  order  accurate  modified  Steger- Warming  flux 
evaluation  method  presented  earlier.  We  performed  perfect  gas  simulations  of  the  flow  at 
several  levels  of  grid  refinement.  The  finest  mesh  used  contained  2048  x  1024  points  in 
the  streamwise  and  wall-normal  directions.  The  results  are  plotted  in  Fig.  10,  which  plots 
heat  transfer  rate  on  the  double-cone  model.  We  see  that  the  heat  transfer  rates  computed 
on  the  1024  x  512  and  2048  x  1024  grids  are  virtually  identical.  Thus,  the  results  obtained 
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Figure  10.  Heat  transfer  rate  for  a  perfect  gas  simulation  computed  on  four 
different  grids;  (b)  shows  an  enlargement  near  the  point  of  separation. 


on  the  1024  x  512  grid  are  converged. 

2.  Importance  of  Time  Integration  Method 

Careful  time  convergence  studies  by  Gaitonde  et  al.  (2002)  showed  that  for  this 
particular  test  case  the  solution  does  not  converge  until  at  least  100  characteristic  flow 
times  have  been  computed.  The  characteristic  flow  time  has  been  defined  as  the  time  it 
takes  for  a  particle  to  travel  the  length  of  the  geometry,  if  it  is  moving  at  the  free-stream 
velocity,  rc har  =  L/uoo-  The  relatively  long  time  scale  of  100xrchar  allows  the  separation 
zone  to  become  fully  established.  In  this  work,  we  simulate  the  flow  for  at  least  150  flow 
times  before  considering  that  the  solution  has  reached  a  steady  state. 

We  achieve  grid  converged  results  on  a  grid  of  1024  x  512  points,  and  the  characteristic 
size  of  the  smallest  cell  on  this  mesh,  Aymin  is  of  the  order  of  0.1  fim.  In  order  to  maintain 
stability  of  the  numerical  scheme,  the  computational  time  allowed  on  such  a  mesh  is  TeXpi  ~ 
CFLexpi  x  Aymin/ttoo,  where  CFLexpi  <  0.4  for  second-order  schemes.  Then,  to  reach  the 
steady  state  solution  with  an  explicit  method  of  time  integration,  Nexp\  —  150Tchar/VeXpi 
iterations  must  be  performed.  The  double-cone  geometry  is  approximately  0.18  m  in 
length,  and  therefore,  Nexp\  is  of  the  order  108.  Performing  such  a  large  number  of  iterations 
would  require  an  enormous  amount  of  physical  time,  thus  making  the  simulation  of  the 
double-cone  flows  impossible.  Therefore,  it  is  necessary  to  employ  an  implicit  method  of 
time  integration,  which  allows  for  a  time  step  much  larger  than  rexpi  to  be  taken. 

3.  Flux  Evaluation  Methods 

Several  flux  evaluation  methods  are  employed  in  this  work,  with  the  purpose  to  ex¬ 
amine  the  effect  of  numerical  dissipation  on  the  solution  of  double-cone  flows.  In  addition 
to  the  modified  Steger- Warming  method  discussed  previously,  the  Roe  method  with  the  H 
entropy  correction  [Sanders  et  al.  (1998)],  the  Lax-Friedrichs  (LF)  method,  the  Harten- 
Lax-van  Leer  (HLL)  and  Harten-Lax-van  Leer  Contact  (HLLC)  methods  are  employed. 
These  methods  are  well  documented  in  the  literature  [Toro  (1997)].  Also,  higher  order  of 
spatial  accuracy  of  the  flux  evaluation  is  achieved  using  the  MUSCL  variable  extrapolation 
method,  where  the  conserved  variables  are  extrapolated  to  the  cell  faces  using  a  second  or¬ 
der  reconstruction  scheme.  Non-linear  limiters  are  commonly  used  to  prevent  oscillations, 
and  the  minmod,  van  Albada,  van  Leer,  and  superbee  limiters  are  employed  in  this  work. 
For  a  discussion  on  slope  limiters  see  Hirsch  (1991).  More  details  regarding  limiters  can 
be  found  in  Laney  (1998). 

The  fluxes  computed  using  the  different  flux  evaluation  methods  must  be  accurately 
linearized  in  time  in  order  to  be  used  in  conjunction  with  an  implicit  solver.  The  lineariza¬ 
tion  procedure  depends  on  the  flux  evaluation  method,  and  in  some  cases  increases  the 
computational  cost  associated  with  the  original  method.  The  linearization  procedure  pre- 
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sented  previously  can  be  applied  directly  to  split  flux  methods,  but  the  same  approach  can 
be  taken  to  obtain  flux  Jacobians  for  other  methods  [Druguet  et  al.  (2003)].  Furthermore, 
the  implicit  equation  remains  unchanged  for  the  inviscid  problem. 

4.  Time  Convergence 

We  are  interested  in  obtaining  the  steady  state  solution  with  each  of  the  different 
numerical  methods,  thus  we  may  start  the  simulation  with  any  physically  acceptable  initial 
condition.  We  typically  start  the  simulation  by  initializing  the  flow  to  the  free-stream 
conditions  everywhere.  However,  we  may  also  start  from  any  solution  already  obtained 
using  any  of  the  flux  evaluation  methods,  and  iterate  until  the  solution  reaches  steady 
state.  This  may  result  in  considerable  savings  in  computational  time. 

Because  we  use  an  implicit  solver,  we  take  the  largest  possible  time  step  in  order 
to  obtain  results  in  the  smallest  computational  time.  However,  if  the  time  step  is  too 
large,  the  computation  does  not  converge  as  fast  to  the  steady  result  as  with  a  smaller 
time  step.  If  the  time  step  is  even  larger,  the  simulation  may  enter  a  limit  cycle,  where  the 
solution  keeps  oscillating  among  states,  and  the  residual  never  decreases.  This  is  commonly 
referred  to  as  ringing  of  the  solution.  This  may  be  particularly  important  for  flows  with 
shock  interactions  and  recirculation  regions.  Therefore,  the  convergence  rate  depends 
highly  on  the  time  step  size.  Also,  simulations  with  different  numerical  flux  evaluation 
methods  converge  differently.  For  example,  the  HLL  and  HLLC  methods  require  more 
iterations  than  the  Roe  scheme  to  converge.  Additionally,  simulations  performed  with  the 
Roe  scheme  and  the  minmod  limiter  converge  much  faster  than  simulations  performed  with 
the  same  scheme  and  less  dissipative  limiters,  such  as  the  van  Albada  and  superbee. 

5.  Results  on  Coarse  Grids 

In  this  subsection  we  present  results  of  simulations  performed  on  the  coarsest  mesh 
(256  x  128)  used  in  this  work.  We  show  the  shock  location  by  plotting  the  T  =  400  K  tem¬ 
perature  contour  in  Fig.  11.  The  numerical  results  were  computed  with  different  second 
order  flux  evaluation  methods  (mod  SW,  LF,  Roe,  HLL,  and  HLLC).  These  computed  re¬ 
sults,  with  the  exception  of  the  one  computed  with  the  modified  Steger- Warming  method, 
use  a  solution  reconstruction  scheme  where  the  slopes  are  limited  with  the  minmod  limiter. 
This  limiter  is  known  to  be  dissipative.  From  this  figure,  we  see  that  the  size  of  the  sep¬ 
aration  zone  depends  on  the  method  used.  The  more  dissipative  schemes  (LF  and  HLL) 
produce  a  much  smaller  separation  zone  than  the  less  dissipative  schemes  (mod  SW  and 
Roe).  The  solution  computed  with  the  modified  Steger- Warming  method  gives  the  largest 
separation  zone  because  the  scheme  used  for  variable  extrapolation  is  much  less  dissipative 
than  the  minmod  limiter.  When  the  Roe  scheme  is  used  with  a  slope  limiter  less  dissipative 
than  minmod,  it  results  in  a  larger  separation  zone.  This  is  shown  in  Fig.  12,  which  plots 
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Figure  11.  Shock  location  (T  =  400  K  contour)  computed  with  different  flux 
evaluation  methods  on  a  256  x  128  grid. 


Figure  12.  Shock  location  (T  =  400  K  contour)  computed  with  the  Roe  scheme 
and  different  slope  limiters  on  a  256  x  128  grid. 
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Figure  13.  Surface  pressure  predictions  on  the  double-cone  computed  with  dif¬ 
ferent  flux  evaluation  methods  on  a  256  x  128  grid. 

the  shock  location  computed  with  different  slope  limiters  for  the  variable  extrapolation. 
Interestingly,  when  we  rank  the  results  from  the  smallest  to  the  largest  size  of  separation, 
the  list  corresponds  to  the  classification  of  the  slope  limiters  according  to  decreasing  dis¬ 
sipation  [Toro  (1997)].  These  results  show  that  the  less  dissipative  the  numerical  scheme, 
the  larger  the  separation  zone.  Furthermore,  these  results  are  reflected  in  the  plots  of  the 
surface  quantities;  the  pressure  is  shown  in  Fig.  13. 

6.  Results  on  Fine  Grids 

In  this  subsection,  we  show  results  obtained  on  a  fine  mesh.  Fig.  14  shows  the  surface 
pressure  computed  on  a  1024  x  512  grid  using  four  different  schemes  (mod.  SW,  Roe  with 
minmod,  HLLC  with  minmod,  and  LF  with  minmod).  These  results  are  much  closer  to 
one  another  than  the  results  obtained  on  the  coarse  grid,  shown  in  Fig.  13.  The  results 
obtained  using  the  Roe  and  HLLC  schemes  with  the  minmod  slope  limiter  are  identical 
on  this  mesh.  The  most  dissipative  scheme  (LF+  minmod)  still  does  not  give  the  same 
results  as  the  other  methods,  because  it  has  not  yet  reached  grid  convergence.  The  size 
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Figure  14.  Surface  pressure  on  the  double-cone  computed  with  different  flux 
evaluation  methods  on  a  1024  x  512  grid. 

of  the  separation  zone  obtained  with  the  less  dissipative  Roe  scheme  is  still  smaller  than 
what  is  obtained  with  the  modified  Steger- Warming  method.  This  is  because  the  minmod 
limiter  is  too  dissipative. 

There  is  noticeable  improvement  in  the  results  when  a  less  dissipative  slope  limiter  is 
used  for  variable  extrapolation,  especially  when  the  solver  itself  is  very  dissipative.  This 
is  shown  in  Fig.  15,  which  plots  the  surface  pressure  computed  using  the  Lax-Friedrichs 
scheme,  and  two  different  limiters  (minmod  and  van  Leer).  Also,  if  we  compare  the  Roe 
and  Lax-Friedrichs  solutions  computed  with  the  van  Leer  limiter  in  both  cases,  we  see  that 
the  results  are  very  similar. 

7.  Evolution  of  the  Separation  Zone  Size 

It  is  very  useful  to  know  how  different  numerical  schemes  perform  by  looking  at  a  sin¬ 
gle  scalar  quantity,  as  opposed  to  distributions  of  surface  quantities.  In  earlier  subsections 
we  showed  how  the  size  of  the  recirculation  zone  depends  greatly  on  the  amount  of  numer¬ 
ical  dissipation  of  the  numerical  scheme.  Thus,  we  can  get  a  measure  of  the  dissipation  by 
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Figure  15.  Surface  pressure  on  the  double-cone  model  computed  with  with  the 
Lax-Friedrichs  scheme  and  different  slope  limiters  on  a  1024  x  512  grid. 

comparing  the  size  of  the  separation  zone  that  is  predicted  when  a  particular  scheme  is 
employed.  Additionally,  we  can  obtain  information  about  grid  convergence  by  plotting 
the  size  of  the  separation  zone  versus  the  size  of  the  mesh.  First,  we  need  to  define 
how  we  measure  physically  the  size  of  the  recirculation  region.  We  define  the  size  of  the 
recirculation  L,  as  the  length  between  the  point  of  detachment  on  the  first  cone,  and  the 
point  of  re-attachment  on  the  second  cone.  It  is  easy  to  obtain  this  information  from 
numerical  solutions  by  plotting  the  surface  skin  friction  coefficient  Cf.  The  detachment 
point  is  where  the  skin  friction  coefficient  becomes  negative,  and  the  re-attachment  point 
is  where  the  Cf  becomes  positive  again.  Note  that  in  certain  cases  there  can  be  multiple, 
small-scale,  bubbles  of  separated  flow  inside  the  large  recirculation  zone. 

Fig.  16  plots  the  size  of  the  recirculation  zone  obtained  with  different  schemes,  versus 
the  square  of  the  inverse  of  the  number  of  points  in  the  streamwise  direction  Nx.  This 
figure  allows  us  to  show  the  evolution  of  the  size  of  the  recirculation  zone  as  the  grid 
is  refined,  when  different  numerical  schemes  are  employed.  The  results  obtained  on  the 
coarsest  grid  are  further  to  the  right.  As  the  grid  is  refined,  1/NX  decreases,  and  the  results 
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Figure  16.  Size  of  the  recirculation  zone  versus  the  square  of  the  inverse  of  the 
number  of  points  in  the  streamwise  direction. 

converge  toward  the  same  value.  When  the  grid  is  refined  to  2048  x  1024  points,  the  result 
obtained  with  the  modified  Steger- Warming  method  remains  the  same  as  on  the  1024  x 
512  grid.  This  shows  that  our  results  are  grid  converged. 

8.  Conclusions 

This  study  shows  that  numerical  predictions  of  hypersonic  flow  over  the  double-cone 
model  are  very  sensitive  to  the  numerical  flux  evaluation  method.  All  schemes  give  the 
same  results  when  grid  convergence  is  reached,  however  on  coarse  meshes  different  flux 
evaluation  schemes  give  different  results.  The  least  dissipative  flux  evaluation  methods 
examined  in  this  work  (modified  Steger- Warming,  Roe  with  van  Leer  limiter,  and  HLL- 
Contact  with  van  Leer  limiter)  give  grid  converged  results  on  a  1024  x  512  grid.  The  more 
dissipative  schemes  (Lax-Friedrichs,  Harten-Lax-van  Leer)  and  schemes  that  employ  the 
minmod  limiter  for  variable  extrapolation  do  not  converge  at  this  level  of  grid  refinement. 

It  is  important  to  note  that  the  grid-converged  results  obtained  on  the  1024  x  512  mesh 


49 


using  the  least  dissipative  method  do  not  match  the  experimental  data.  Most  importantly, 
the  size  of  the  separation  zone  appears  to  be  over-predicted  by  the  most  numerically 
accurate  solution.  This  is  not  due  to  numerical  error,  and  it  is  because  of  inadequate 
modeling  of  the  flow  physics,  as  we  will  see  in  the  following  section. 

Because  the  double-cone  flow  takes  a  long  time  to  reach  its  steady  state,  large  amounts 
of  physical  time  must  be  simulated.  Therefore,  an  implicit  time  integration  method  is 
required  for  simulating  these  flows  within  a  reasonable  time  frame. 

Finally,  the  dissipation  properties  of  a  numerical  scheme  can  be  assessed  and  compared 
with  other  methods  by  using  the  double-cone  flow.  In  particular,  by  measuring  the  size  of 
the  recirculation  zone  from  a  given  numerical  solution,  a  direct  assessment  can  be  made 
of  the  quality  of  the  numerics  employed  in  the  simulation.  Furthermore,  this  can  be  done 
without  having  to  perform  a  complete  grid  convergence  study. 
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Characterization  of  Experiments 

In  the  fourth  section  we  presented  results  of  simulations  over  two  double-cone  flows, 
and  we  made  comparisons  with  the  experimental  data.  We  observed  good  overall  agreement 
between  numerical  predictions  and  experimental  measurements,  with  the  exception  of  the 
heat  transfer  rate  on  the  first  cone  near  the  tip.  This  discrepancy  was  common  to  both  cases 
presented,  and  these  results  are  consistent  with  simulations  performed  by  other  researchers. 
Therefore,  the  purpose  of  this  section  is  to  present  a  study  of  the  source  of  this  discrepancy. 
We  do  this  by  examining  the  way  the  experiments  are  performed,  and  in  particular  how 
the  free-stream  conditions  are  obtained.  We  also  examine  possible  non-continuum  effects 
on  solid  boundaries. 

We  assess  the  accuracy  of  the  nominal  test-section  conditions  by  using  CFD  to  com¬ 
pute  the  flow  in  the  LENS  facility  nozzle,  including  the  effects  of  vibrational  nonequilib¬ 
rium.  We  also  investigate  the  accuracy  of  the  conventional  no-slip  boundary  condition  at 
solid  walls  for  these  relatively  low-density  flows.  Finally,  we  use  our  computed  test-section 
conditions  to  study  the  effect  of  slight  free-stream  non-uniformity  on  the  computed  sepa¬ 
ration  zone  size.  When  all  these  effects  are  included  in  the  simulation  of  the  double-cone 
flow  of  interest,  improved  agreement  between  computations  and  experiments  is  obtained. 

1.  Vibrational  Nonequilibrium 

The  specification  of  the  free-stream  conditions  is  an  important  aspect  of  hypersonic 
flow  experiments.  In  a  reflected  shock  hypersonic  wind  tunnel  such  as  LENS  Leg  I,  high 
temperature  and  pressure  gas  is  expanded  from  the  reservoir  conditions  through  a  con¬ 
toured  nozzle  to  high  Mach  number.  It  is  well  known  that  the  vibrational  modes  of  the  gas 
may  freeze  near  the  nozzle  throat  conditions.  Vibrational  freezing  is  particularly  impor¬ 
tant  in  nitrogen  at  relatively  low  pressure  reservoir  conditions,  as  is  the  case  in  the  present 
experiments. 

In  this  subsection,  we  discuss  how  vibrational  freezing  affects  the  inferred  test  con¬ 
ditions  in  the  LENS  facility  for  this  particular  class  of  experiments.  In  the  experiments, 
calibration  runs  are  performed  before  the  tests  to  verify  the  uniformity  of  the  flow  sec¬ 
tion.  Subsequently,  the  actual  experiments  are  performed  at  the  same  nominal  conditions 
as  the  calibration  runs.  In  an  impulse  facility  such  as  LENS,  the  pressure  and  enthalpy 
in  the  reservoir  region  behind  the  reflected  shock  are  known.  The  Pitot  pressure  in  the 
test-section  can  be  measured  easily,  as  well  as  the  heat  transfer  rate  to  a  reference  probe. 
Additional  Pitot  pressure  measurements  are  made  during  the  experiment  to  verify  that 
the  test  conditions  are  consistent  with  the  calibration  run. 

To  determine  the  free-stream  conditions,  a  quasi-one-dimensional  code  is  run  using 
the  measured  reservoir  stagnation  conditions.  Although  the  exact  geometrical  specification 
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of  the  nozzle  is  known,  the  effective  area  ratio  is  not  known  due  to  the  wall  boundary  layer 
displacement.  Therefore,  the  code  is  run  to  the  point  where  the  computed  Pitot  pressure 
matches  the  pressure  measured  by  the  Pitot  probe  at  the  exit-plane  of  the  nozzle.  Then, 
the  test  conditions  are  taken  to  be  those  given  by  the  quasi-one-dimensional  code.  This 
analysis  assumes  vibrational  equilibrium  during  the  expansion,  and  has  been  successfully 
used  for  a  long  time,  but  usually  in  air  at  pressures  higher  than  considered  in  the  present 
experiments. 


To  illustrate  how  the  free-stream  conditions  are  inferred,  consider  the  Rayleigh  Pitot 
pressure  formula 
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which  in  the  hypersonic  limit  reduces  to 


_  7  +  1  f  (7  +  l)2 
2  V  47 


7  +  1  ( (7  +  l)2 
2q  y  4q 


P00M 


Poouoo 


where  p0 2  is  the  measured  Pitot  pressure  and  00  denotes  the  free-stream  conditions.  Thus, 
the  Pitot  pressure  is  effectively  a  measurement  of  the  kinetic  energy  per  unit  volume  in 
the  free-stream. 


Now  consider  an  ideal,  adiabatic  expansion  from  the  reservoir  to  hypersonic  conditions. 
The  total  enthalpy  of  the  reservoir  is  converted  to  free-stream  kinetic  energy  conserved 
according  to 

ha  —  CpToo  4-  2Uoo  —  2U°° 

where  we  have  assumed  that  \u^  »  CpT^.  Therefore,  given  the  reservoir  enthalpy,  we 
have  the  value  of  u ^  in  the  inviscid  portion  of  the  nozzle  flow.  Then,  with  a  measurement 
of  the  Pitot  pressure,  we  effectively  measure  the  free-stream  density  for  the  ideal  expansion. 


Now  consider  what  happens  if  there  is  vibrational  energy  frozen  in  the  flow  during  the 
expansion.  We  have 


where  e*  is  the  vibrational  energy  per  unit  mass  frozen  in  the  flow.  Thus,  the  kinetic 
energy  will  be  lower  than  in  an  ideal  expansion,  and  to  obtain  the  same  measured  Pitot 
pressure,  the  free-stream  density  must  be  larger.  Thus,  the  effect  of  vibrational  freezing  in 
the  nozzle  is  to  lower  the  axial  velocity  and  increase  the  density  relative  to  an  equilibrium 
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analysis.  To  first-order,  the  surface  pressure  scales  with  p0 oit^,  and  the  heat  transfer  rate 
with  PoqU^,.  Thus,  we  do  not  expect  vibrational  freezing  to  affect  the  measured  surface 
pressure,  but  it  should  reduce  the  measured  heat  transfer  rate  due  to  the  lower  value  of 

^OO* 

In  this  analysis  we  have  assumed  that  the  Rayleigh  Pitot  pressure  formula  is  valid 
under  conditions  of  vibrational  nonequilibrium.  We  have  simulated  the  flow  over  the  Pitot 
probe  used  in  the  experiments  using  a  finite-rate  relaxation  model  for  nitrogen.  These 
calculations  show  that  under  the  conditions  of  interest,  the  Rayleigh  Pitot  pressure  formula 
gives  the  correct  value  of  Pitot  pressure,  and  the  above  analysis  is  valid. 

2.  Nozzle  Flow  Simulations 

The  preceding  elementary  analysis  shows  that  vibrational  freezing  reduces  the  mea¬ 
sured  heat  transfer  rate.  However,  a  more  complete  analysis  is  required  to  determine  how 
much  vibrational  freezing  occurs  and  to  assess  the  additional  non-ideal  effects  in  the  nozzle 
flow.  Thus,  we  have  performed  a  series  of  CFD  simulations  of  the  nozzle  flows.  In  these 
simulations,  we  have  used  the  numerical  method  presented  in  the  third  section  with  the 
Baldwin-Lomax  (1978)  and  the  Spalart-Allmaras  (1992)  turbulence  models  for  the  nozzle 
wall  boundary  layer. 

The  nozzle  grids  are  constructed  directly  from  the  CAD  files  of  the  nozzle  contours. 
Typical  grids  use  1788  points  in  the  axial  direction  and  128  points  in  the  surface  normal 
direction,  with  stretching  at  the  surface.  The  nozzle  throats  are  elongated,  making  the 
specification  of  the  sonic  line  difficult.  Thus,  the  flow  is  computed  from  the  constant 
diameter  driven  tube  section,  through  the  throat,  and  to  the  test-section.  Thermo-chemical 
equilibrium  was  assumed  at  the  inflow,  and  subsonic  inflow  conditions  were  implemented 
so  that  total  enthalpy  is  conserved  in  the  computed  flow. 

Initial  simulations  of  the  nozzle  flow  showed  a  significant  over-prediction  of  the  center- 
line  Pitot  pressure.  This  is  caused  by  the  Baldwin-Lomax  turbulence  model  over-predicting 
the  boundary  layer  displacement  thickness  at  the  nozzle  exit-plane.  This  results  in  higher 
velocity  in  the  inviscid  region,  and  also  high  Pitot  pressure.  The  strong  favorable  pres¬ 
sure  gradient  in  the  nozzle  reduces  the  effect  of  turbulent  transport  through  the  boundary 
layer,  reducing  the  rate  of  boundary  layer  growth.  This  mechanism  is  not  included  in  any 
algebraic  turbulence  models.  Thus,  we  recognized  that  the  results  obtained  using  the  al¬ 
gebraic  turbulence  model  were  not  realistic.  Therefore,  we  implemented  the  one-equation 
turbulence  model  of  Spalart  and  Allmaras,  with  the  compressibility  correction  suggested 
by  Catris  and  Aupoix  (2000),  to  perform  additional  simulations.  This  model  allows  for 
suppression  and  amplification  of  turbulence  based  on  the  local  flow  characteristics.  Im¬ 
plementation  of  this  model  requires  the  solution  of  an  additional  conservation  equation. 
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Figure  17.  Pitot  pressure  in  the  nozzle  exit-plane  compared  with  nozzle  compu¬ 
tation  for  Run  1302  (calibration  for  Run  35). 

We  simulated  the  flow  in  the  nozzle  using  this  model  and  found  that  the  simulation  matched 
the  measured  Pitot  pressure  profile.  This  is  shown  in  Fig.  17,  which  plots  the  measured 
Pitot  pressure  for  the  calibration  run  for  Run  35  (p0  =  3.55  MPa,  h0  =  3.71MJ/kg)  with 
the  results  of  two  nozzle  simulations. 

We  ran  a  nozzle  simulation  with  the  actual  Run  35  reservoir  conditions  (jp0  —  3.77MPa, 
ha  =  3.83MJ/kg).  From  the  simulation,  we  obtain  the  free-stream  conditions  at  the  test- 
section.  The  results  of  this  simulation  are  tabulated  in  Table  2  along  with  those  provided 
by  Holden  (2000)  using  the  equilibrium  quasi-one-dimensional  analysis.  Note  that  the 
vibrational  temperature  is  predicted  to  freeze  at  2562  K,  which  is  close  to  the  throat  tem¬ 
perature.  The  resulting  dynamic  pressure  and  kinetic  energy  flux  are  also  given.  The 
computed  value  of  the  dynamic  pressure  is  lower  than  that  derived  from  the  nominal  con¬ 
ditions  because  we  have  matched  the  Pitot  pressure  profile,  rather  than  a  single  mean  value 
of  Pitot  pressure.  Also  note  that  PoouL  12.4%  lower  in  the  nonequilibrium  simulations 
due  to  the  6.7%  lower  value  of  Pitot  pressure  and  the  6.2%  lower  value  of  axial  velocity. 
A  similar  analysis  was  carried  out  for  Run  28,  and  the  results  are  given  in  Table  3.  In  this 
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case,  the  predicted  Pitot  pressure  is  2.5%  higher  and  the  resulting  kinetic  energy  flux  is 
reduced  by  only  2.6%. 


Case 

Nominal 

Nonequilibrium 

(K) 

138.9 

98.27 

TVo o  (K) 

138.9 

2562 

Poo  (g/m3) 

0.5515 

0.5848 

Uoo  (m/s) 

2713 

2545 

Moo 

11.30 

12.59 

Pcoulc  (Pa)  4059  3788 

Pooi&>  (W/cm2)  1101  964 


Table.  2  Nominal  and  computed  test-section  conditions  for  Run  35  reservoir  con¬ 
ditions. 


Case  Nominal  Nonequilibrium 


Too  (K) 

185.6 

140.0 

Tv oo  (K) 

185.6 

2589 

Poo  (g/m3) 

0.6545 

0.7372 

Uoo  (m/s) 

2664 

2538 

Moo 

9.50 

10.5 

Po ou^o  (Pa) 

4645 

4749 

Pooulo  (W/cm2) 

1237 

1205 

Table.  3  Nominal  and  computed  test-section  conditions  for  Run  28  reservoir  con¬ 
ditions. 

3.  Results 

In  the  previous  sections,  we  quantified  the  effect  of  nonequilibrium  in  the  nozzle  and 
obtained  nonequilibrium  free-stream  conditions  by  performing  simulations  of  the  nozzle 
flow.  In  this  subsection,  we  use  the  nonequilibrium  conditions  as  free-stream  specification 
to  recompute  the  flows  over  double-cone  models  and  make  comparisons  with  the  measured 
surface  quantities.  We  show  that  under  nonequilibrium  free-stream  conditions  the  heat 
transfer  rate  to  the  first  cone  is  reduced.  However,  this  reduction  does  not  result  in  com¬ 
plete  agreement  with  the  experimental  measurements.  We  also  consider  non-continuum 
effects,  and  in  particular  the  effect  of  partial  accommodation  of  incident  particles  at  solid 
boundaries.  We  examine  the  effect  of  velocity  slip  and  temperature  “jump”  at  the  surface  of 
the  model,  as  well  as  possible  jump  of  vibrational  energy.  We  employ  the  standard  Maxwell 
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a) 


Figure  18.  Heat  transfer  rate  and  surface  pressure  for  nominal  Run  35  conditions. 
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model  for  velocity  slip  and  temperature  jump,  with  a  simple  extension  to  include  jump  of 
vibrational  energy  as  described  previously.  In  this  model,  a,  ot  are  the  accommodation 
coefficients  for  velocity  and  temperature  respectively,  and  are  taken  to  be  0.85.  The  value 
of  the  vibrational  energy  accommodation  coefficient  is  not  well  known,  and  a  value  of  0.5 
was  used  initially.  This  value  of  ov  results  in  nearly  complete  vibrational  accommodation. 
However,  Meolans  et  at  (2002)  suggest  adiabatic  conditions  (av  —  0),  and  Black  et 
at  (1974)  measure  values  of  the  order  10-3  for  stainless  steel.  We  show  that  when 
the  slip  model  is  included  in  the  simulations  under  nominal  free-stream  conditions  the 
results  are  very  similar  to  the  no-slip  results.  Under  nonequilibrium  conditions  and  with 
<jv  =  1.0  x  10-3  we  see  a  substantial  decrease  in  heat  transfer  rate  to  the  first  cone.  Finally, 
we  investigate  how  small  levels  of  test-section  flow  non-uniformity  in  the  free-stream  affect 
the  flow.  We  focus  on  Run  35,  and  present  more  limited  results  for  Rim  28. 

The  effect  of  slip  at  Run  35  nominal  conditions  is  shown  in  Fig.  18.  In  this  com¬ 
putation,  the  accommodation  coefficients  for  velocity  and  temperature  were  0.85,  while  a 
value  of  0.5  was  used  for  vibrational  energy.  We  see  that  slip  has  little  effect,  with  the 
largest  difference  being  the  heat  transfer  to  the  first  cone,  which  decreases  by  less  than 
3%.  Therefore,  slip  is  not  by  itself  responsible  for  the  over-prediction  of  heat  transfer  to 
the  first  cone. 

We  observe  a  more  significant  decrease  in  heat  transfer  rate  at  the  first  cone  when  the 
computed  nonequilibrium  test-section  conditions  are  employed,  as  shown  in  Fig.  19.  The 
heat  transfer  rate  decreases  by  9.8%  on  the  first  cone,  and  the  surface  pressure  decreases  by 
4.7%  in  that  region,  better  matching  the  experimental  data.  The  pressure  in  the  separation 
zone  is  also  reduced  (by  6.3%),  improving  the  agreement  with  the  experiment.  Note  that 
the  separation  point  and  the  locations  of  the  peak  heat  transfer  and  pressure  do  not  change 
substantially.  However,  the  heat  transfer  rate  and  pressure  on  the  last  2  cm  of  the  second 
cone  are  reduced  by  11%  and  8.1%,  respectively.  We  also  employed  the  slip  model  using  the 
nonequilibrium  test-section  conditions  with  the  same  accommodation  coefficients.  There 
is  no  significant  difference  in  surface  properties,  with  heat  transfer  being  reduced  by  less 
than  2.0%.  Thus,  including  both  vibrational  nonequilibrium  in  the  nozzle  and  surface  slip, 
the  heat  transfer  to  the  first  cone  is  reduced  by  12%,  and  the  difference  is  still  greater  than 
the  experimental  uncertainty. 

Note  that  at  the  predicted  nonequilibrium  test-section  conditions,  the  vibrational 
temperature  is  highly  elevated.  Therefore,  the  gradient  of  vibrational  energy  at  the  room- 
temperature  surface  is  large,  and  the  modeling  of  vibrational  energy  accommodation  is 
important.  Theory  and  experiments  suggest  that  the  vibrational  energy  accommodation 
coefficient  should  be  much  smaller  than  the  value  of  0.5  used  in  the  previous  simulations. 
Therefore,  we  simulated  Run  35  with  the  nonequilibrium  conditions  and  av  =  0.001;  the 
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Figure  19.  Heat  transfer  rate  and  surface  pressure  for  nonequilibrium  Run  35 
conditions. 


Figure  20.  Heat  transfer  rate  for  nonequilibrium  Run  35  conditions  and  av  = 
0.001. 

results  are  plotted  in  Fig.  20.  The  pressure  distribution  and  separation  zone  size  are 
unaffected.  But  the  predicted  heat  transfer  rate  on  the  first  cone  decreases  by  an  additional 
7.0%,  to  within  the  5%  uncertainty  in  the  data.  Interestingly,  the  predicted  heat  transfer 
to  the  back  half  of  the  second  cone  is  now  below  the  measurements. 

Clearly  vibrational  nonequilibrium  plays  an  important  role  in  these  flows.  It  affects 
the  test-section  conditions  by  reducing  the  flux  of  kinetic  energy.  And  because  of  weak 
vibrational  accommodation,  a  large  fraction  of  the  free-stream  vibrational  energy  is  not 
transferred  to  the  surface  near  the  tip. 

An  additional  non-ideal  effect  is  the  weak  non-uniformity  in  the  test-section  flow. 
Fig.  21  plots  the  computed  Mach  number  contours  in  the  test-section,  with  the  double¬ 
cone  flow  field  superimposed  at  the  test  location.  To  assess  the  effect  of  the  flow  non¬ 
uniformity,  we  extract  flow  conditions  along  a  line  that  corresponds  to  the  location  of  the 
inflow  boundary  condition  for  the  double-cone  simulation.  Fig.  22  compares  this  calculation 
with  the  previous  av  —  0.001  result.  The  free-stream  non-uniformity  reduces  the  size  of  the 
separation  zone  by  3.6  mm,  shifting  the  heat  transfer  peak  upstream  in  better  agreement 
with  the  data.  Also  note  improved  agreement  on  either  side  of  the  peak. 
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Figure  21.  Computed  Mach  number  contours  in  the  test-section  and  with  double 
cone  flow  solution  superimposed. 


x  (cm) 

Figure  22.  Heat  transfer  rate  for  non-uniform  Run  35  conditions  and  av  =  0.001 
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Now  consider  the  analysis  of  Run  28.  At  nominal  test  conditions,  slip  has  little 
effect,  as  in  Run  35.  Fig.  23  shows  that  when  we  use  the  computed  nonequilibrium  test- 
section  conditions,  the  heat  transfer  rate  to  the  first  cone  is  reduced  by  just  2.1%,  which  is 
consistent  with  the  small  (2.6%)  reduction  in  the  kinetic  energy  flux  for  this  case.  Fig.  23 
also  shows  that  the  nonequilibrium  conditions  produce  a  significantly  larger  separation 
zone,  apparently  worsening  the  agreement  with  experiment.  Fig.  24  shows  that  vibrational 
energy  slip  decreases  the  heat  transfer  rate  by  11%,  bringing  it  within  about  8.7%  of  the 
measurements.  However,  the  size  of  the  separation  zone  increases.  Fig.  25  shows  the  effect 
of  test-section  non-uniformity  for  Run  28.  The  simulations  now  agree  well  with  the  data, 
due  to  a  smaller  predicted  separation  zone  and  a  slightly  reduced  heat  transfer  rate  to  the 
first  cone. 

These  results  show  that  including  vibrational  energy  freezing  in  the  nozzle,  weak 
vibrational  accommodation  to  the  model  surface,  and  test-section  non-uniformity  gives 
excellent  agreement  between  computations  and  experiments.  It  should  be  noted  that 
these  experiments  were  run  in  pure  nitrogen  at  low  pressures,  which  tends  to  enhance 
the  effects  of  vibrational  nonequilibrium.  Also  the  non-uniformity  of  the  test-section  is 
likely  increased  by  the  low  pressure  conditions,  which  result  in  thicker  boundary  layers 
and  non-ideal  nozzle  performance.  Therefore,  it  is  likely  that  under  the  usual  operating 
conditions  of  the  LENS  facility  (air  at  orders  of  magnitude  larger  pressure)  these  effects 
will  be  diminished.  However  further  studies  are  required  to  quantify  this  statement. 

4.  Under-resolved  Solutions 

In  this  section  we  discuss  the  effect  of  grid  resolution  on  the  double-cone  flows,  and 
show  that  if  the  flow  field  is  under-resolved  it  is  possible  to  obtain  good  agreement  with 
experiments  for  the  wrong  reason.  The  computations  presented  in  the  previous  section 
are  computationally  demanding  because  the  recirculation  zone  develops  slowly  and  a  large 
physical  time  must  be  simulated  to  reach  steady  state.  In  addition,  large  finely-spaced 
grids  must  be  used  to  ensure  grid-converged  results.  For  example,  to  capture  the  initial 
boundary  layer  growth  at  the  cone  tip,  our  first  grid  point  is  located  0.1  yum  from  the 
surface. 

Let  us  study  how  the  solution  varies  when  the  flow  is  under-resolved.  Fig.  26  plots  the 
heat  transfer  rate  at  nominal  and  nonequilibrium  conditions  for  Run  28  using  a  256  x  128 
grid.  The  coarse  grid  solutions  produce  smaller  separation  zones,  and  as  a  result  the 
nonequilibrium  calculation  appears  to  give  the  “correct”  extent  of  separation.  However, 
Fig.  27  shows  that  when  we  include  the  additional  effects  of  small  vibrational  energy 
accommodation  and  flow  non-uniformity,  the  coarse  grid  no  longer  agrees  so  well  with  the 
experiments.  This  is  an  excellent  illustration  of  the  importance  of  performing  a  rigorous 
grid  convergence  study,  preferably  without  knowledge  of  the  data. 
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Figure  23.  Heat  transfer  rate  and  surface  pressure  for  nonequilibrium  Run  28 
conditions. 
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Figure  24.  Heat  transfer  rate  for  nonequilibrium  Run  28  conditions  and  av  = 
0.001. 


Figure  25.  Heat  transfer  rate  for  non-uniform  Run  28  conditions  and  av  =  0.001. 
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Figure  26.  Heat  transfer  rate  for  (a)  nominal  and  (b)  nonequilibrium  Run  28 
conditions  computed  on  the  coarse  grid  (256  x  128). 


Figure  27.  Heat  transfer  rate  for  non-uniform  Run  28  conditions  with  crv  =  0.001 
computed  on  the  coarse  grid  (256  x  128). 

5.  Results  for  Low  Density  Experiments 

In  this  section  we  present  simulation  results  of  double-cone  experiments  performed 
at  very  low  density  nitrogen  and  specific  enthalpy  approximately  equal  to  that  of  Run 
35.  These  low  density  double-cone  experiments  were  performed  following  the  initial  blind 
comparisons  with  experimental  data,  and  were  designed  to  be  used  as  validation  cases  for 
particle-based  methods  such  as  the  Direct  Simulation  Monte  Carlo  (DSMC)  method.  The 
reason  for  performing  additional  double-cone  experiments  with  substantially  lower  free- 
stream  densities  is  to  make  the  problem  more  tractable  for  the  DSMC  codes  by  reducing 
the  computational  cost.  The  cost  of  particle-based  simulations  increases  with  the  number 
of  computational  particles  that  must  be  simulated.  The  existing  double-cone  experiments 
were  designed  to  be  in  the  continuum  regime.  This,  combined  with  the  large  amount  of 
physical  time  that  must  be  simulated  for  the  double-cone  flow  to  reach  steady  state  make 
the  DSMC  simulation  prohibitively  expensive. 

We  present  simulation  results  for  two  low  density  cases:  Run  4  and  Run  5.  The 
experimental  conditions  for  these  cases  were  designed  to  result  in  low  density  free-streams, 
with  approximately  1/2  and  1/3  the  free-stream  density  of  Run  35.  These  are  Run  4  and 
Run  5,  corresponding  to  the  1/2  density  and  1/3  density  conditions  respectively.  The  free- 
stream  conditions  for  these  cases  are  listed  in  Table  4.  The  free-stream  conditions  were 
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obtained  assuming  vibrational  nonequilibrium  in  the  nozzle  flow,  and  we  see  that  there  is 
vibrational  energy  frozen  in  the  flow  as  was  the  case  for  Run  28  and  Run  35. 


Case 

Run  4 

“1/2  density” 

Run  5 

“1/3  density” 

Moo 

15.90 

15.27 

Rem  (103/m) 

109.0 

55.5 

Too  (K) 

46.58 

52.99 

Tv  oo  (K) 

2027.1 

2128.4 

Twall  (K) 

294.48 

292.78 

Poo  (10-3  kg/m3) 

0.3020 

0.1606 

Uoo  (m/s) 

2213.6 

2266.4 

Table.  4  Free-stream  and  wall  conditions  for  the  “1/2  density”  (Run  4)  and  “1/3 
density”  (Run  5)  25°-55°  sharp  double-cone  LENS  experiments. 

Fig.  28  shows  the  results  of  simulations  for  Run  4,  and  the  heat  transfer  rate  to 
the  body  is  plotted  for  two  simulations  along  with  the  experimental  measurements.  The 
simulation  under  free-stream  conditions  with  vibrational  nonequilibrium  included  shows 
good  agreement  with  the  measurements  everywhere,  with  the  exception  of  the  heat  transfer 
rate  to  the  first  cone.  This  result  is  consistent  with  our  previous  analysis  of  Run  35  and 
Run  28.  When  the  effects  of  slip  are  included  in  the  simulation  and  the  accommodation 
coefficient  for  vibrational  energy  is  set  to  a  very  small  value  (0.001),  the  heat  transfer  rate 
to  the  first  cone  is  reduced  by  approximately  10%,  and  good  agreement  is  observed  for  the 
1/2  density  case.  It  is  also  interesting  to  note  that  there  is  very  good  agreement  between 
both  simulations  and  the  measurements  on  either  side  of  the  peak  in  heat  transfer  rate 
and  on  the  second  cone.  However,  it  is  surprising  that  the  location  of  the  peak  heating  is 
predicted  correctly,  while  the  separation  zone  size  is  under-predicted  by  about  2  mm.  This 
is  probably  due  to  slight  non-uniformity  in  the  test-section,  as  was  the  case  for  Run  35. 

Fig.  29  plots  similar  heat  transfer  rate  results  for  Run  5.  The  separation  zone  size  is 
predicted  correctly,  and  the  separation  region  is  much  smaller  for  this  case  than  the  higher 
density  case  (Run  4).  This  is  because  the  free-stream  Reynolds  number  is  much  lower  for 
this  case,  and  this  result  is  consistent  with  results  of  previous  studies.  The  heat  transfer 
rate  to  the  first  cone  is  over-predicted  by  about  10%  with  the  no-slip  simulation.  However, 
the  heat  transfer  rate  to  the  first  cone  is  lowered  and  is  within  the  experimental  5% 
uncertainty  when  slip  effects  with  weak  vibrational  energy  accommodation  to  the  surface 
are  modeled. 

From  these  simulation  results  we  see  that  the  effect  of  vibrational  nonequilibrium  in 
the  free-stream  and  non-continuum  effects  on  the  surface  are  present  in  double-cone  flows 
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Figure  28.  Heat  transfer  rate  for  nonequilibrium  Run  4  conditions  with  av  =  0.5 
(solid  line)  and  crv  =  0.001  (dashed  line). 


Figure  29.  Heat  transfer  rate  for  nonequilibrium  Run  5  conditions  with  crv  =  0.5 
(solid  line)  and  crv  =  0.001  (dashed  line). 
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of  different  free-stream  Mach  and  Reynolds  numbers.  This  shows  that  we  have  correctly 
identified  the  source  of  the  initial  discrepancies  and  our  analysis  is  valid. 

6.  Conclusions 

In  order  to  understand  why  discrepancies  were  observed  between  initial  numerical 
simulations  and  the  experiments  presented  previously,  we  investigated  effects  of  vibrational 
nonequilibrium,  slip  at  the  model  surface,  and  weak  flow  non-uniformity.  We  simulated 
the  flow  within  the  facility  nozzle,  and  found  that  under  the  low  pressure  conditions  of 
the  experiment,  the  vibrational  energy  freezes  near  the  throat  temperature.  This  reduces 
the  kinetic  energy  flux  in  the  nozzle  test-section,  reducing  the  heat  transfer  to  the  model. 
However,  the  heat  transfer  rate  was  still  over-predicted.  Therefore,  we  considered  the 
failure  of  the  no-slip  boundary  condition  at  the  model  surface.  We  found  that  the  computed 
heat  transfer  rate  is  sensitive  to  the  vibrational  energy  accommodation  coefficient.  Using 
experimental  values,  the  heat  transfer  to  the  first  cone  is  predicted  to  within  the  5% 
measurement  uncertainty.  Finally,  we  used  the  results  of  our  nozzle  flow  field  predictions 
as  inflow  to  the  double-cone  flow  solution.  These  simulations  accounted  for  small  levels 
of  flow  non-uniformity  in  the  test-section,  and  show  an  improvement  in  the  size  of  the 
separation  zone  and  the  heat  transfer  rate  to  the  second  cone. 

In  this  study,  we  also  illustrated  the  importance  of  conducting  careful  grid  convergence 
studies.  Because  the  computed  separation  zone  size  depends  on  the  grid  resolution,  it  is 
possible  to  obtain  spurious  “agreement”  with  experiments  due  to  poor  grid  resolution  and 
inadequate  modeling  of  the  flow  physics. 
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High  Enthalpy  Double-Cone  Experiments 

In  this  section  we  present  simulations  of  double-cone  flow  experiments  at  high  enthalpy 
free-stream  conditions  performed  at  the  LENS  facility.  These  experiments  were  conducted 
through  a  joint  computational  and  experimental  effort,  based  on  the  methodology  for 
characterizing  experiments  that  was  presented  in  the  previous  section.  This  approach  was 
successful  in  developing  improved  test  conditions  and  physical  models  for  low  enthalpy 
nitrogen  experiments. 

The  present  experiments  are  of  high  enthalpy  hypersonic  flows,  characteristic  of  flight 
at  a  Mach  number  of  about  14.  New  experiments  have  been  performed  using  the  LENS 
Leg  I  tunnel  at  these  conditions  for  both  nitrogen  and  air.  The  flow  in  the  facility  nozzle  is 
simulated,  and  flows  over  double-cones  are  computed  at  the  nominal  test  conditions,  and 
at  the  conditions  predicted  by  the  nozzle  simulations. 

1.  Design  of  High  Enthalpy  Experiments 

In  previous  work  [Nompelis  et  al.  (2003a)],  we  performed  simulations  of  the  flow  of 
air  in  the  nozzle  near  the  maximum  LENS  Leg  I  enthalpy  conditions  ( hQ  =  11.5  MJ/kg, 
p0  —  59  MPa).  We  systematically  lowered  the  reservoir  pressure,  keeping  the  total  en¬ 
thalpy  constant,  and  simulated  the  flow  in  the  nozzle  from  the  driven  section  to  the  nozzle 
exit-plane.  We  obtained  computed  free-stream  conditions  for  four  different  reservoir  condi¬ 
tions;  the  existing  conditions  of  p0  —  7380  psi  and  three  additional  cases  (6000,  4000,  and 
2000  psi).  Then,  we  simulated  the  flow  over  the  double-cone  model  under  the  correspond¬ 
ing  computed  test-section  conditions  that  we  obtained  from  the  nozzle  flow  simulations.  In 
order  to  determine  whether  flows  under  these  proposed  conditions  are  laminar,  we  exam¬ 
ined  the  shear  layer  that  forms  at  the  edge  of  the  separation,  which  is  the  most  unstable 
region  of  the  flow.  We  compared  the  shear  layer  Reynolds  numbers  with  the  Birch  and 
Keyes  criterion  for  transition  in  free  shear  layers  formed  by  shock  interactions  [Birch  and 
Keyes,  (1972)].  Based  on  this  analysis,  we  found  that  the  double-cone  flows  generated 
under  the  two  lowest  stagnation  pressure  conditions  will  likely  be  laminar.  These  results 
are  consistent  with  previous  laminar  double-cone  flows  under  low  enthalpy  conditions. 

Experiments  were  performed  at  the  proposed  lower  pressure  conditions  in  the  LENS  I 
facility  and  measurements  were  made  over  the  re-instrumented  double-cone  model.  Under 
these  low-pressure  conditions,  we  anticipated  that  there  would  be  some  level  of  chemical 
reaction  and  vibrational  nonequilibrium  in  the  free-stream.  Thus,  the  effects  of  chemical 
reactions  and  energy  relaxation  in  the  nozzle  were  taken  into  consideration  when  inferring 
the  free-stream  conditions.  Determination  of  the  free-stream  conditions  involves  perform¬ 
ing  measurements  of  the  Pitot  pressure  at  the  nozzle  exit-plane  using  Pitot  probes.  We 
performed  simulations  including  real-gas  effects  of  the  flow  around  the  Pitot  probes  and 
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found  that  the  Pitot  pressure  is  measured  accurately  by  the  probe  under  high  enthalpy 
conditions.  Once  the  Pitot  pressure  at  the  test-section  is  known,  a  quasi  one-dimensional 
code  is  used  to  expand  the  reservoir  conditions  and  match  the  Pitot  pressure.  Finite-rate 
chemistry  with  vibrational  relaxation  of  the  gas  has  been  integrated  in  this  algorithm. 
The  vibrational  state  of  the  gas  is  computed  based  on  a  vibrational  state-specific  model 
for  nitrogen  and  a  vibrational  equilibrium  model  for  air.  This  assumption  of  vibrational 
equilibrium  is  based  on  previous  experience  at  Calspan  with  high  enthalpy  shock  tunnel 
flows  and  because  of  the  difficulty  of  predicting  coupled  vibrational  relaxation  and  recom¬ 
bination  in  expanding  air.  The  free-stream  conditions  for  high  enthalpy  experiments  in 
pure  nitrogen  and  air  are  listed  in  Table  5.  Note  that  there  is  a  substantial  degree  of 
vibrational  nonequilibrium  in  the  nitrogen  flow.  Also  in  the  air  flow,  there  is  some  NO  and 
O  in  the  free-stream  as  predicted  by  this  analysis,  and  therefore  chemical  energy  is  stored 
in  the  gas. 


Case 

Run42 

(nitrogen) 

Run  46 
(nitrogen) 

Run  43 
(air) 

Moo 

11.52 

11.54 

8.87 

Rem  (103/m) 

333.4 

439.4 

306.6 

Too  (K) 

268.7 

281.7 

576.0 

Tv  oo  (K) 

3160 

3072 

576.0 

Twall  (K) 

294.7 

296.3 

296.2 

Poo  (g/m3) 

1.468 

1.958 

2.134 

Uoo  (m/s) 

3849.3 

3946.9 

4218.1 

cn2 

1.00000 

0.99842 

0.73704 

co2 

- 

- 

0.17160 

cno 

- 

- 

0.06477 

Cn 

0.00000 

0.00158 

0.00000 

Co 

- 

- 

0.02659 

Table.  5  Nominal  free-stream  and  wall  conditions  for  the  25-55°  sharp  double-cone 
high  enthalpy  experiments. 

2.  High  Enthalpy  Double-Cone  Flows 

We  simulated  the  flow  over  the  double-cone  model  under  the  nominal  conditions  for 
the  nitrogen  and  air  cases.  A  schematic  of  the  high  enthalpy  double-cone  flow  in  ni¬ 
trogen  (Run  42)  is  shown  in  Fig.  30.  The  interacting  shocks  create  a  massive  region 
of  separation  that  extends  to  about  half  the  length  of  the  first  cone.  The  shear  layer 
that  forms  at  the  edge  of  the  separation  is  about  6  cm  in  length.  A  supersonic  jet  is 
formed  near  the  surface  of  the  second  cone.  Interestingly,  there  are  two  regions  behind  the 
bow  shock  where  the  flow  is  subsonic.  Fig.  31  plots  the  surface  pressure  and  heat  transfer 
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Figure  30.  Schematic  of  double-cone  flow-field  at  high  enthalpy  nitrogen  free- 
stream  conditions  (Run  42). 

rate  for  Run  42.  In  general,  agreement  between  the  simulation  and  the  experimental  mea¬ 
surements  is  very  good.  The  size  of  the  separation  zone  is  almost  exactly  reproduced,  as 
well  as  the  location  of  the  pressure  peak  on  the  second  cone.  Surprisingly,  the  pressure 
is  under-predicted  along  the  first  cone  and  in  the  separation  region,  by  10%  and  5%  re¬ 
spectively.  The  heat  transfer  is  also  under-predicted  along  the  first  cone  and  at  the  peak. 
Also,  the  heat  transfer  peak  in  the  experiments  appears  to  be  wider  than  the  simulation 
predicts.  The  effect  of  slip  on  the  heat  transfer  rate  is  more  pronounced  along  the  first 
cone,  while  the  pressure  is  unaffected.  This  result  is  consistent  with  the  low  enthalpy  Run 
35  results. 

Fig.  32  plots  the  results  for  Run  46,  which  is  nitrogen  at  a  40%  higher  stagnation 
pressure  than  Run  42.  The  higher  Reynolds  number  results  in  a  slightly  larger  separation 
zone,  as  expected.  We  see  that  the  agreement  is  good.  However  we  under-predict  the  size 
of  the  separation  zone,  and  in  general  the  predicted  heat  transfer  and  pressure  is  slightly 
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Figure  31.  Surface  pressure  (a)  and  heat  transfer  rate  to  the  model  (b)  for  high 
enthalpy  nitrogen  free-stream  conditions  (Run  42). 


low.  It  is  unlikely  that  this  difference  is  due  to  transition  to  turbulence,  because  the  flow 
at  the  forebody  is  certainly  laminar,  and  heat  transfer  rate  is  predicted  to  be  low  there. 

Fig.  33  shows  a  comparison  between  experimental  measurements  and  simulation  re¬ 
sults  for  the  air  case  (Run  43).  The  size  of  the  separation  zone  is  not  predicted  correctly, 
although  the  basic  flow  has  been  captured  by  the  CFD.  The  heat  transfer  rate  is  under¬ 
predicted  by  about  30%  in  the  shock  interaction  region,  and  the  heat  transfer  rate  to  the 
second  cone  is  also  significantly  under-predicted  (by  20-30%).  In  general,  the  size  of  the 
separation  zone  affects  the  location  of  the  peak  and  the  strength  of  the  interaction.  In 
this  case,  the  small  separation  zone  results  in  having  the  peak  in  heat  transfer  slightly 
upstream  from  what  the  experiments  show.  Surface  catalysis  is  not  responsible  for  this 
large  disagreement  as  shown  by  the  fully  catalytic  simulation  (catalytic  efficiency  was  set 
to  unity  at  the  wall). 

The  flowfields  generated  under  the  new  experimental  conditions  appear  to  be  fully 
laminar.  We  computed  the  shear  layer  Reynolds  numbers  based  on  the  Birch  and  Keyes 
definition,  and  the  results  are  plotted  in  Fig.  34.  In  this  figure,  the  points  represent  the 
Reynolds  number  at  different  locations  along  the  shear  layer,  computed  based  on  the  state 
of  the  gas  on  the  high-speed  side:  Res  =  pvs/p,  where  s  is  the  length  along  the  shear 
layer.  This  dimensionless  quantity  remains  under  20,000  for  Run  42  and  under  10,000  for 
Run  43,  and  this  is  well  under  the  Birch  and  Keyes  criterion  for  transition.  The  shear  layer 
Reynolds  numbers  for  the  fully  laminar  low  enthalpy  Run  35  case  are  plotted  for  reference. 

From  these  results  we  see  that  the  high  enthalpy  nitrogen  and  air  double-cone  experi¬ 
ments  are  fully  laminar  and  the  free-stream  conditions  are  calculated  very  accurately.  Also, 
the  results  of  the  numerical  simulation  are  mostly  within  the  5%  experimental  uncertainty 
for  this  effectively  blind  comparison.  In  view  of  the  discrepancies  observed  initially  for  the 
low  enthalpy  experiments  presented  previously,  this  is  encouraging,  particularly  given  the 
uncertainty  introduced  by  the  air  chemistry  models  employed.  It  should  be  noted  that  the 
magnitude  of  the  peak  in  heat  transfer  rate  is  not  reproduced,  and  the  difference  is  sig¬ 
nificant.  This  discrepancy  is  common  among  all  high  enthalpy  flow  simulations  presented 
here. 

3.  Nozzle  Flow  Simulations 

In  this  subsection  we  present  nozzle  flow  simulations  performed  under  the  experimental 
conditions  of  Run  42  and  Run  43.  Our  objective  is  to  provide  additional  characterization 
of  the  free-stream  at  the  high  enthalpy  LENS  conditions.  The  nozzle  flow  is  simulated  from 
the  driven  section,  through  the  nozzle  throat  and  into  the  contoured  nozzle.  The  reservoir 
conditions  are  taken  from  direct  measurements  in  the  driven  section,  where  thermochemi¬ 
cal  equilibrium  is  assumed.  We  use  a  541  x  148  grid  to  represent  the  nozzle,  and  the  grid 
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Figure  33.  Surface  pressure  (a)  and  heat  transfer  rate  (b)  to  the  model  for  high 
enthalpy  air  free-stream  conditions  (Run  43). 
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Figure  34.  Shear  layer  Reynolds  numbers  at  the  edge  of  the  separation  region 
calculated  based  on  the  Birch  and  Keyes  definition  for  Run  42  and  Run  43.  The 
low  enthalpy  Run  35  result  is  plotted  for  reference. 

is  constructed  from  the  nozzle  geometry  specification.  In  previous  studies  of  the  flow 
inside  the  LENS  facility  nozzle  at  low  enthalpy  conditions,  we  used  the  Baldwin-Lomax 
turbulence  model  for  the  nozzle  wall  boundary  layer.  This  turbulence  model  did  not 
reproduce  the  boundary  layer  thickness  observed  in  the  experiments.  However,  when  we 
employed  the  Spalart-Allmaras  model  for  wall  turbulence  the  simulation  matched  the  Pitot 
pressure  profile  at  the  nozzle  exit-plane.  In  the  present  studies,  we  also  employ  the  Spalart- 
Allmaras  model  for  wall  turbulence  [Spalart  and  Allmaras,  (1992)]  with  the  compressibility 
correction  from  Catris  and  Aupoix  (2000). 

Fig.  35  compares  the  computed  exit-plane  Pitot  pressure  profile  with  the  probe  rake 
measurements  for  the  calibration  case  of  Run  42  (Run  187).  The  numerical  predictions 
match  the  measured  values  very  well,  without  any  adjustable  constants  in  the  calcula¬ 
tion.  Fig.  36  plots  the  velocity  profiles  at  different  locations  along  the  nozzle  wall  plotted 
in  wall  variables.  We  see  that  the  turbulence  model  closely  resembles  the  logarithmic 
law  of  the  wall  in  van  Driest  coordinates  in  the  region  of  the  throat  (the  subsonic  region 
and  the  sonic  point),  as  well  as  the  expansion  region  (M  =  2).  This  is  a  verification 
of  the  turbulence  model  implementation  in  the  CFD.  The  profile  progressively  deviates 
from  the  log-law  as  the  flow  expands  to  the  test-section.  The  more  advanced  turbulence 
model  captures  the  reduction  in  the  boundary  layer  thickness  due  to  the  favorable  pressure 
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Figure  35.  Predicted  and  measured  Pitot  pressure  profiles  for  the  calibration  run 
in  nitrogen  (Run  187). 


Figure  36.  Velocity  profiles  of  the  nozzle  wall  boundary  layer  from  the  Run  42 
calibration  case  plotted  in  wall  variables. 
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Figure  37.  Heat  transfer  rate  predictions  and  measurements  at  the  wall  for  Run 
187  and  Run  188. 


Figure  38.  Predicted  and  measured  Pitot  pressure  profiles  for  the  calibration  run 
in  air  (Run  188). 
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gradient.  Fig.  37  plots  the  heat  transfer  to  the  nozzle  wall  during  the  expansion  under 
the  calibration  run  conditions.  The  numerical  results  match  the  measured  values,  and 
the  trend  is  predicted  correctly,  indicating  that  the  boundary  layer  is  turbulent  in  the 
region  where  measurements  could  be  made  in  the  nozzle.  This  is  contrasted  to  results  of 
previous  simulations  in  which  the  algebraic  model  for  wall  turbulence  was  used.  Thus, 
the  use  of  the  Spalart-Allmaras  turbulence  model  makes  the  nozzle  simulations  a  more 
predictive  tool.  The  Pitot  pressure  profile  for  the  air  calibration  run  (Run  188)  is  plotted 
in  Fig.  38  along  with  the  probe  measurements.  The  baseline  (vibrational  equilibrium  and 
nominal  reaction  rates)  simulation  does  not  predict  the  profile  correctly,  and  is  within 
10%  of  the  measured  values.  However,  when  we  force  the  vibrational  energy  to  be  in 
equilibrium  during  the  expansion,  agreement  with  measurements  is  improved.  In  a  similar 
fashion,  if  we  force  the  chemical  reaction  rates  to  be  faster  by  an  order  of  magnitude  for 
the  vibrational  equilibrium  simulations,  we  observe  a  15%  increase  in  the  Pitot  pressure 
at  the  centerline.  Thus,  uncertainties  in  vibrational  relaxation  and  reaction  rate  affect  the 
predicted  test-section  conditions.  The  computed  test-section  conditions  obtained  with  the 
baseline  simulations  are  tabulated  in  Table  6. 


Case 

Run  42 
(nitrogen) 

Run  43 
(air) 

Too  (K) 

268.1 

338.6 

Tv  oo  (K) 

3021.8 

2103.8 

Poo  (10~3  kg/m3) 

1.486 

1.832 

u0 o  (m/s) 

3886.5 

4126.4 

cn2 

0.99765 

0.74422 

co2 

- 

0.15099 

cno 

- 

0.04898 

Cn 

0.00235 

0.00000 

Co 

- 

0.05580 

Table.  6  Computed  test-section  conditions  for  the  25-55°  sharp  double-cone  high 
enthalpy  experiments. 

4.  Simulations  at  Computed  Conditions 

In  this  subsection,  we  present  simulations  of  the  double-cone  flows  using  the  computed 
conditions  presented  in  the  previous  subsection.  Fig.  38  shows  the  computed  pressure  and 
heat  transfer  rate  to  the  model  under  the  computed  conditions  of  Run  42.  We  see  that 
the  pressure  and  heat  transfer  rate  to  the  first  cone  have  increased  by  about  2%,  but 
the  overall  agreement  has  not  improved.  This  is  not  surprising,  given  that  the  computed 
test-section  conditions  are  not  very  different  from  the  nominal  conditions  listed  in  Table 

5.  This  is  because  vibrational  nonequilibrium  was  included  in  the  quasi-one-dimensional 
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Figure  38.  Predicted  surface  pressure  (a)  and  heat  transfer  rate  (b)  for  Run  42 
computed  test-section  conditions. 
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Figure  39.  Predicted  surface  pressure  (a)  and  heat  transfer  rate  (b)  for  Run  43 
computed  test-section  conditions. 
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calculation  used  to  obtain  the  free-stream  conditions  for  the  nitrogen  case. 

Results  of  simulations  of  the  air  case  (Run  43)  under  computed  test-section  conditions 
are  plotted  in  Fig.  39.  A  fully  catalytic  surface  model  was  employed  for  these  simulations. 
The  pressure  and  heat  transfer  rate  are  reduced  by  about  10-15%  under  conditions  for 
vibrational  and  thermal  equilibrium.  These  predictions  substantially  deviate  from  the 
measured  values,  and  also  the  separation  zone  is  under-predicted.  At  computed  conditions 
obtained  with  the  fast  reaction  rates  for  the  nozzle  flow,  the  numerical  predictions  are  very 
close  to  those  at  the  nominal  conditions. 

5.  Conclusions 

We  performed  detailed  numerical  simulations  of  high  enthalpy  nitrogen  and  air  exper¬ 
iments  from  the  LENS  facility.  These  experiments  were  performed  guided  by  the  results 
of  numerical  simulation  and  evaluation  of  previous  high  enthalpy  experiments.  Also,  the 
nominal  experimental  free-stream  conditions  provided  by  the  experimentalists  included 
the  effect  of  vibrational  nonequilibrium.  After  the  experiments  were  completed,  we  re¬ 
evaluated  the  free-stream  conditions  by  performing  numerical  simulations  of  the  nozzle 
flow  under  the  exact  reservoir  conditions  of  the  experiments.  We  employed  the  Spalart- 
Allmaras  turbulence  model  for  all  simulations  of  the  nozzle  flow.  We  were  able  to  verify 
our  implementation  of  the  turbulence  model  by  making  comparisons  with  the  heat  trans¬ 
fer  measurements  to  the  nozzle  wall  during  the  calibration  runs  in  both  air  and  nitrogen. 
Numerical  results  were  in  good  agreement  with  the  Pitot  pressure  measurements  at  the 
nozzle  exit-plane  for  the  calibration  runs.  With  this  turbulence  model,  we  did  not  have 
to  adjust  the  simulation  in  order  to  match  the  Pitot  pressure  profile  at  the  exit-plane. 
Therefore,  the  methodology  that  we  developed  based  on  the  low  enthalpy  experiments  for 
free-stream  characterization  can  be  used  for  experiments  at  high  enthalpy. 

We  simulated  flows  over  the  double-cone  model  at  high  enthalpy  conditions.  We 
compared  numerical  predictions  under  nominal  and  computed  test-section  conditions  with 
the  experimental  data.  The  simulations  are  in  agreement  with  the  measurements  for  the 
pure  nitrogen  cases,  however  the  peak  in  heat  transfer  rate  and  the  width  of  the  peak  are 
significantly  under-predicted  at  higher  Reynolds  numbers.  This  is  indicative  of  inaccurate 
modeling,  either  in  the  determination  of  the  free-stream  conditions,  or  the  chemistry  models 
employed  in  the  simulations.  We  did  not  obtain  as  good  agreement  for  the  air  case.  We 
performed  simulations  where  we  increased  the  rate  of  reaction  both  in  the  nozzle  flow  and 
in  the  flow  over  the  model.  In  this  way,  we  were  able  to  reproduce  the  Pitot  pressure 
profile  in  the  test-section  when  vibrational  equilibrium  was  enforced  in  the  simulation. 
Interestingly,  simulations  under  the  corresponding  test-section  conditions  produced  results 
very  similar  to  those  of  the  baseline  simulation.  We  were  unable  to  improve  agreement 
with  measurements  for  these  high  enthalpy  experiments. 
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Design  of  Experiments  in  the  Expansion  Tunnel 

In  the  previous  section  we  presented  simulations  of  high  enthalpy  air  and  nitrogen 
experiments  performed  using  the  LENS  Leg  I  tunnel.  These  experiments  were  performed 
at  low  pressures,  such  that  the  resulting  flow-fields  were  laminar  and  steady.  The  results 
showed  good  overall  agreement  for  nitrogen  flows,  but  significant  differences  between  mea¬ 
surements  and  predictions  were  observed  for  air  flows.  This  is  not  a  surprising  result. 
The  pressures  at  which  the  experiments  were  performed  were  very  low,  in  order  to  ensure 
that  the  flow  over  the  model  is  laminar.  At  such  low  pressures,  the  free-stream  generated 
by  the  reflected  shock  tunnel  is  in  vibrational  nonequilibrium  and  is  weakly  dissociated. 
Therefore,  some  degree  of  uncertainty  is  introduced  because  of  this  effect.  Therefore,  an 
alternative  method  of  generating  high  Mach  number  free-streams  at  low  pressure  has  been 
considered,  and  this  involves  using  an  expansion  tunnel  instead  of  a  reflected  shock  tunnel. 

A  new  hypervelocity  tunnel  has  been  developed  at  the  LENS  facility,  which  utilizes  an 
expansion  tube  to  generate  high  Mach  number  free-streams,  and  is  designated  as  LENS- 
X.  In  this  subsection  ,  we  present  a  preliminary  study  of  a  high  enthalpy  experiment 
performed  using  the  LENS-X  expansion  tunnel.  Computational  fluid  dynamics  is  used  to 
simulate  the  entire  experiment.  This  includes  simulations  of  the  flow  inside  the  expansion 
tube  and  nozzle  after  the  primary  diaphragm  breaks,  and  also  the  flow  over  the  double¬ 
cone  model.  First,  we  provide  a  short  description  of  the  facility,  and  discuss  results  of  the 
one-dimensional  analysis  performed  for  two  cases.  We  proceed  by  presenting  axisymmetric 
simulations  of  the  flow  inside  the  expansion  tunnel  assuming  laminar  and  fully  turbulent 
flow  for  one  of  the  conditions.  We  obtain  free-stream  conditions  by  performing  simulations 
of  the  nozzle  flow,  and  finally,  we  simulate  the  flow  over  the  model  using  these  conditions. 

1.  Expansion  Tube  Experiments 

Previous  high  enthalpy  experiments  at  the  LENS  facility  were  performed  using  the 
LENS  Leg  I  tunnel,  operated  in  reflected  shock  mode.  In  these  experiments,  the  test  gas 
in  the  reservoir  region  behind  the  diaphragm  is  shock  heated  twice  and  stagnates  before 
it  accelerates  into  the  nozzle.  This  allows  the  gas  to  reach  its  equilibrium  composition  at 
the  stagnation  temperature  in  the  reservoir  region.  As  the  flow  accelerates  through  the 
throat  and  into  the  nozzle,  the  gas  cools  rapidly  during  the  expansion  before  it  reaches 
the  test-section.  Depending  on  the  flow  conditions,  this  process  may  result  in  a  partially 
dissociated  or  nonequilibrium  free-stream.  This  is  likely  to  be  the  case  for  the  low  pressure 
experiments  of  interest.  Therefore,  an  alternative  method  of  obtaining  realistic  high  en¬ 
thalpy  free-streams  using  the  existing  experimental  facilities  was  devised  by  making  some 
modifications  to  the  LENS  facility.  A  section  was  added  to  the  LENS  Leg  II  tunnel,  which 
allows  it  to  be  operated  as  an  expansion  tunnel.  The  details  can  be  found  in  Holden  (2004). 
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Figure  40.  Schematic  of  an  expansion  tube  and  its  x  —  t  diagram. 

The  added  section  makes  a  transition  from  the  expansion  tube  to  the  existing  divergent 
nozzle.  This  type  of  arrangement  forces  the  gas  to  undergo  a  steady  expansion  inside  the 
divergent  nozzle,  following  the  unsteady  expansion  it  experiences  in  the  expansion  tube 
[Lukasiewicz  (1973)]. 

The  expansion  tube  is  a  device  by  which  high  enthalpy  cold  flows  can  be  obtained  ex¬ 
perimentally.  The  device  consists  of  a  long  tube,  separated  by  a  primary  and  a  secondary 
diaphragm  as  shown  in  Fig.  40.  The  three  compartments  contain  the  driver  and  driven 
gas  at  different  temperatures  and  pressures.  The  driver  gas  is  at  high  pressure  and  tem¬ 
perature,  and  is  in  the  unperturbed  region  1  as  indicated  on  the  x  —  t  diagram  of  Fig.  40. 
The  driven  gas  occupies  two  compartments,  the  second  of  which  is  at  a  substantially  lower 
pressure  than  the  first;  these  are  regions  2  and  3  in  Fig.  40.  When  the  primary  diaphragm 
breaks  at  the  desired  pressure,  the  driver  gas  accelerates  into  the  second  compartment, 
compressing  the  driven  gas,  and  sending  a  shock  wave  traveling  down  the  tube,  as  seen  in 
the  figure.  Region  4  is  a  region  of  uniform  state  behind  the  primary  shock,  and  is  the  test 
gas.  When  the  primary  shock  breaks  the  second  diaphragm,  it  sends  a  fast-moving  shock 
down  the  expansion  region.  The  total  enthalpy  is  conserved  and  thus  the  gas  behind  the 
shock  in  region  5  is  at  very  high  temperature.  The  gas  in  region  6  has  been  accelerated  to 
high  velocity,  but  is  at  a  relatively  low  temperature.  Therefore,  if  air  is  placed  in  region  2 
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as  the  test  gas,  it  can  be  set  into  motion  at  high  enthalpy  and  maintained  at  a  non-reacted 
state  (region  6)  using  this  experimental  setup. 

The  expansion  tunnel  facility  LENS-X  consists  of  a  nearly  48  m  long,  0.6  m  (24 
in)  diameter  tube  with  a  nozzle  attached  to  it.  Experiments  have  been  performed  in  the 
LENS-X  facility  for  two  initial  conditions.  The  driver  gas  for  these  experiments  was  helium 
and  the  test  gas  was  air.  The  three  compartments  were  filled  according  to  the  conditions 
tabulated  in  Table  7.  In  the  experiments,  the  primary  diaphragm  is  actively  broken,  and 
the  primary  shock  heats  the  test  gas  in  the  583  inch  long  driven  section,  The  second 
diaphragm  breaks  when  the  shock  strikes  it,  and  the  test  gas  is  accelerated  behind  the 
contact  surface  that  is  trailing  the  newly  formed  shock  wave. 


Compartment 

Case 

T(°R) 

A 

V  (psi) 

Case 

T(°R) 

B 

V  (psi) 

Driver  gas 

770.0 

1200 

630.0 

2000 

Driven  gas 

529.3 

3.000 

527.9 

1.470 

Expansion  section 

529.3 

0.0480 

527.9 

0.0059 

Table.  7  Initial  experimental  conditions  for  LENS-X. 

In  the  following  subsections  we  present  results  from  a  one-dimensional  simulation 

of  this  process  for  both  experimental  conditions.  We  use  computational  fluid  dynamics 

to  perform  this  analysis,  assuming  an  inviscid  flow  with  finite  rate  chemical  reactions 

and  vibrational  relaxation.  The  conditions  of  case  B  result  in  vibrational  nonequilibrium 

£ 

inside  the  expansion  tube.  This  is  not  the  case  for  conditions  A,  where  the  test  gas  is  in 
equilibrium.  Subsequently,  we  perform  two-dimensional  laminar  and  turbulent  simulations 
of  the  flow  in  the  expansion  tube  for  case  A.  From  these  simulations  we  extract  free-stream 
conditions  for  the  flow  over  the  model. 

2.  Computational  Procedure 

Simulations  of  the  expansion  tube  are  performed  using  the  numerical  method  pre¬ 
sented  in  section  three.  We  solve  the  conservation  equations  for  a  six-species  helium-air 
mixture.  An  additional  equation  is  solved  along  with  the  conservation  equations  for  the 
unsteady  RANS  simulations  of  the  expansion  tube  and  the  nozzle  flows.  The  Spalart- 
Allmaras  one-equation  turbulence  model  is  employed  in  its  original  form,  and  the  equation 
is  solved  in  conservation  law  form.  The  MUSCL  scheme  for  variable  extrapolation  with  a 
minmod  limiter  is  also  used  selectively  in  this  work.  We  state  explicitly  where  this  scheme 
is  employed. 
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3.  One-Dimensional  Analysis 

In  this  subsection  we  present  results  from  the  one-dimensional  simulation  of  the  flow 
inside  the  expansion  tube.  We  show  that  this  type  of  analysis  is  very  useful  for  assessing 
the  state  and  composition  of  the  test  gas  during  the  experiment.  For  these  simulations, 
a  first-order  accurate  stencil  was  used  on  a  10,000  point  mesh  representing  the  48  m  long 
tube. 

Consider  the  experimental  setup  based  on  the  initial  conditions  of  case  A  mentioned 
previously.  The  primary  diaphragm  breaks  and  we  allow  the  flow  to  evolve  by  integrating 
the  equations  in  a  time-accurate  manner.  This  process  is  illustrated  in  Fig.  41,  which 
plots  the  velocity  and  temperatures  on  the  left  vertical  axis,  and  the  mass-fractions  of  the 
species  on  the  right  axis,  versus  axial  length  measured  from  the  primary  diaphragm.  The 
gas  behind  the  primary  shock  is  heated  and  is  set  in  motion.  This  is  seen  in  Fig.  41a  as 
the  sharp  increase  in  translational-rotational  temperature  and  velocity.  When  the  shock 
ruptures  the  second  diaphragm,  a  fast-moving  shock  forms  in  the  expansion  region,  and  the 
test  gas  goes  through  a  rapid  expansion.  The  temperature  behind  the  shock  wave  reaches 
levels  at  which  chemical  reactions  take  place.  This  is  indicated  by  the  increase  in  the 
mass-fractions  of  N,  O,  and  NO  and  the  corresponding  decrease  in  N2  and  O2  (Fig.  41b). 
The  effect  of  numerical  dissipation  is  visible  by  the  contact  surfaces,  which  are  no  longer 
sharp.  This  is  not  critical  for  the  purpose  of  this  analysis.  There  are  two  interesting  things 
to  note  for  this  simulation.  First,  the  small  region  of  moving  test  gas  behind  the  shock  is  at 
a  uniform  state  -  this  is  region  6  in  Fig.  40.  Secondly,  for  this  particular  initial  condition, 
the  test  gas  is  in  its  initial  composition  and  in  equilibrium  before  it  enters  the  nozzle.  The 
hot,  partially  reacted  gases  behind  the  shock  will  enter  the  nozzle  and  pass  over  the  model 
geometry.  Therefore,  it  is  critical  for  the  measurements  over  the  model  surface  to  be  taken 
when  the  test  gas  is  flowing  over  the  geometry. 

We  performed  the  same  analysis  for  case  B.  The  simulation  proceeds  in  the  same 
manner,  and  the  flow  is  very  similar  to  that  of  case  A.  The  test  gas  experiences  an  expansion 
behind  the  fast-moving  shock,  while  the  air  in  the  expansion  section  dissociates  in  the 
high  temperature  region.  Interestingly,  the  test  gas  is  not  in  equilibrium  as  it  leaves 
the  expansion  section  and  enters  the  nozzle  (Fig.  42).  Based  on  this  result,  we  did  not 
investigate  case  B  any  further,  and  focused  our  work  on  case  A. 

4.  Two-Dimensioned  Simulations 

From  the  one-dimensional  analysis  we  were  able  to  calculate  the  approximate  time- 
varying  state  of  the  gas  as  it  enters  the  nozzle.  However,  it  is  possible  that  the  flow 
exhibits  strong  two-dimensionality,  both  in  the  expansion  tube  and  in  the  nozzle  itself. 
Thus,  a  more  detailed  two-dimensional  simulation  was  required  in  order  to  have  an  accurate 
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b) 


Figure  41.  Evolution  of  the  1-D  flow  in  the  expansion  tube  after  rupture  of  (a) 
the  primary  diaphragm  and  (b)  the  secondary  diaphragm  for  case  A. 
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Figure  42.  Snapshot  of  the  1-D  flow  in  the  expansion  tube  for  case  B. 

time-varying  input  specification  for  the  nozzle  flow.  We  carefully  constructed  a  8000  x  160 
point  mesh  to  represent  the  entire  tube,  with  very  close  stretching  near  the  wall.  After 
having  performed  preliminary  axisymmetric  simulations  on  coarser  meshes,  we  concluded 
that  this  size  mesh  would  produce  reasonably  accurate  results  without  making  the  simu¬ 
lations  extremely  expensive.  In  this  subsection,  we  show  that  the  flow  inside  the  tube  is 
most  likely  turbulent,  and  that  with  a  laminar  simulation  we  are  unable  to  get  realistic 
results. 

We  performed  simulations  assuming  a  laminar  axisymmetric  flow  inside  the  tube.  For 
this  simulation,  the  higher-order  reconstruction  for  the  evaluation  of  the  inviscid  fluxes 
was  based  on  the  MUSCL  variable  extrapolation  scheme,  and  the  minmod  limiter  was 
used.  We  used  this  type  of  reconstruction  because  with  this  method,  extrapolation  is 
performed  individually  for  each  of  the  conserved  variables.  This  is  appropriate  for  the 
unsteady  turbulent  simulations,  where  the  turbulent  viscosity  may  exhibit  gradients  in 
regions  where  other  flow  quantities  are  uniform  and  vice  versa.  Thus,  the  formulation  is 
consistent  between  laminar  and  turbulent  simulations.  The  flow  behind  the  primary  shock 
is  mainly  uniform,  but  exhibits  small  oscillations  and  vortical  structures  near  the  wall. 
The  presence  of  these  oscillations  is  probably  the  cause  of  the  weak  two-dimensionality 
that  we  observed  in  these  flows. 
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Figure  43.  Snapshot  of  the  laminar  (top)  and  fully  turbulent  (bottom)  flows 
inside  the  expansion  tube  -  temperature  contours  are  shown. 

The  Reynolds  number  near  the  wall  region  of  the  flow  behind  the  primary  diaphragm 
is  in  the  order  of  a  few  million  per  meter.  Thus,  the  boundary  layer  that  forms  near  the 
wall  is  likely  to  be  turbulent.  Therefore,  we  performed  an  unsteady  RANS  simulation  of 
the  nozzle  flow  assuming  that  the  flow  is  fully  turbulent.  The  turbulent  simulation  predicts 
a  relatively  uniform  flow  profile  in  the  radial  direction,  throughout  the  entire  length  of  the 
tube.  Small  oscillations  were  observed  shortly  after  the  primary  diaphragm  rupture,  but 
they  did  not  persist  after  the  second  diaphragm  was  ruptured. 

Fig.  43  shows  contours  of  translational  temperature  for  the  laminar  and  turbulent 
simulations,  just  before  the  primary  shock  enters  the  nozzle.  The  test  gas  is  traveling  to 
the  right,  and  is  situated  to  the  left  of  the  hot  gas.  We  can  see  that  there  is  some  two- 
dimensionality  in  the  laminar  simulation,  while  the  turbulent  simulation  resembles  “plug 
flow.”  Fig.  44  plots  flow  variables  at  the  centerline  of  the  tube  for  the  snapshots  shown  in 
Fig.  43.  At  the  centerline  the  flows  are  very  similar,  with  the  exception  of  the  presence  of 
small  oscillations  in  the  laminar  simulation.  In  the  turbulent  simulation,  there  is  a  region 
of  the  flow  behind  the  primary  shock,  where  the  gas  is  at  a  uniform  state.  This  is  not  the 
case  for  the  laminar  simulation,  where  the  variation  of  the  temperature  and  velocity  in 
that  region  are  significant.  This  is  attributed  to  the  presence  of  small  disturbances  near 
the  advancing  front  during  the  evolution  of  the  flow  in  the  tube.  In  view  of  these  results, 
we  concluded  that  a  laminar  simulation  gives  substantially  different  results  than  the  more 
realistic  turbulent  one.  Thus,  it  is  not  possible  to  obtain  a  realistic  time-varying  input 
specification  for  the  nozzle  flow  using  a  laminar  simulation  of  the  expansion  tube. 

5.  Nozzle  Flow  Simulations 

In  this  subsection  we  present  results  from  the  simulation  of  the  nozzle  flow.  The 
flow  inside  the  nozzle  is  assumed  to  be  fully  turbulent,  and  is  computed  in  a  way  sim¬ 
ilar  to  the  flow  inside  the  expansion  tube.  We  extended  the  computational  domain 
ahead  of  the  primary  shock  shown  in  Fig.  43,  by  attaching  a  2000  x  160  point  mesh 
representing  the  nozzle  geometry.  In  this  way,  we  were  able  to  simply  continue  the 
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a) 


b) 


Figure  44.  State  variables  and  species  mass-fractions  at  the  centerline  of  the  tube 
for  laminar  (left)  and  fully  turbulent  (right)  flow  snapshots  of  Fig.  43. 
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Figure  45.  Computed  and  measured  time  history  of  the  Pitot  pressure  (left)  and 
computed  wall  pressure  (right)  near  the  exit  of  the  expansion  section. 

evolution  of  the  flow,  and  thus  we  effectively  provided  a  time  dependent  input  specification 
to  the  nozzle  flow.  This  simulation  was  performed  using  the  MUSCL  scheme. 

In  this  simulation,  we  allowed  the  flow  to  evolve  until  the  driver  gas  reached  the  test- 
section.  During  the  experiments,  measurements  are  made  of  the  Pitot  pressure  and  wall 
pressure  at  the  exit  of  the  expansion  tube.  These  measurements  can  be  compared  directly 
with  the  computed  Pitot  and  static  pressures,  which  are  shown  in  Fig.  45.  The  Pitot  and 
wall  pressures  are  plotted  versus  the  time  measured  from  when  the  primary  shock  reaches 
the  nozzle  section.  We  see  the  sharp  increase  in  Pitot  pressure  when  the  primary  shock 
passes,  followed  by  another  increase  when  the  test  gas  arrives  (Fig.  45a).  We  also  see  the 
increase  in  the  wall  pressure  when  the  test  gas  passes  (Fig.  45b). 

The  test  gas  expands  to  a  high  Mach  number  at  the  centerline  just  upstream  of  the 
exit-plane.  Then,  the  gas  recompresses  to  a  Mach  number  of  about  7  as  it  enters  the 
test-section.  We  should  emphasize  that  when  the  test  gas  enters  the  nozzle,  it  undergoes 
a  steady  expansion.  This  steady  expansion  is  desired  because  it  allows  for  a  steady  nozzle 
flow  to  be  established.  For  this  particular  experimental  condition,  we  predict  the  flow 
inside  the  nozzle  to  remain  steady  for  an  excess  of  3  msec. 

Fig.  46  plots  the  Pitot  pressure  at  the  exit-plane  of  the  nozzle  versus  the  time  from 
when  the  primary  shock  enters  the  nozzle.  The  primary  shock  takes  approxi  mately  4 
msec  to  reach  the  exit-plane  of  the  nozzle.  The  oscillations  that  follow  are  due  to  the 
structure  of  the  advancing  front  that  consists  of  the  hot  gas.  When  the  test  gas  reaches 
the  exit-plane,  the  Pitot  pressure  at  the  centerline  remains  constant  at  about  55  psi, 
as  indicated  by  the  plateau  in  Fig.  46.  A  direct  estimate  of  the  conditional  test  time  can 
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Figure  46.  Computed  time  history  of  the  Pitot  pressure  at  the  nozzle  exit-plane. 

be  made  from  this  plateau,  which  measures  approximately  3  msec.  Note  that  this  Pitot 
measurement  is  at  the  centerline.  Therefore  an  assessment  must  be  made  of  the  uniformity 
of  the  flow  in  the  test-section  by  examining  the  radial  distribution  of  flow  quantities. 
Furthermore,  the  distribution  of  flow  quantities  in  the  radial  direction  must  be  examined 
over  the  duration  of  the  useful  test  time.  We  do  this  by  extracting  flow  profiles  in  the 
radial  direction  for  three  instances  marked  by  the  arrows  on  Fig.  46,  and  indicated  by 
ti  —  5.7  msec,  1 2  =  7.2  msec,  and  <3  =  8.7  msec.  The  radial  distributions  of  axial  velocity, 
Pitot  pressure,  and  Mach  number  are  plotted  in  Fig.  47.  The  profiles  show  that  the  flow 
is  relatively  uniform  inside  a  circular  region  of  about  0.5  m  (20  in)  in  diameter.  This  is 
because  of  the  large  size  of  the  facility.  Also,  we  see  that  the  time-varying  profiles  start 
by  exhibiting  substantial  non-uniformity  immediately  following  the  arrival  of  the  test  gas. 
This  is  because  of  the  structure  of  the  advancing  front. 

With  this  study  of  the  flow  we  were  able  to  get  an  estimate  of  the  test  time  that  can 
be  achieved  using  this  setup,  and  also  calculate  the  state  of  the  gas  at  the  test-section. 
The  computed  free-stream  conditions  are  tabulated  in  Table  8. 

6.  Double-Cone  Flow 

We  simulated  the  flow  over  the  25°-55°  double-cone  model  with  free-stream  conditions 
as  given  in  Table  8.  We  used  a  mesh  of  512  x  256  points  in  the  streamwise  and  wall-normal 
directions.  The  numerical  method  employed  was  the  same  as  what  we  used  to  simulate  the 
low  and  high  enthalpy  double-cone  flows  presented  previously.  This  particular  numerical 
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method  along  with  a  mesh  of  this  size  is  expected  to  give  results  accurate  to  about  1%  in 
terms  of  the  size  of  the  separation  zone.  The  computed  flowfield  exhibited  a  very  large 
region  of  separation.  The  existence  of  such  a  large  region  of  separated  flow  is  an  indication 
that  this  flow  might  be  transitional.  This  is  probably  because  the  free-stream  density  is  an 
order  of  magnitude  larger  than  that  of  the  earlier  LENS  Leg  I  double-cone  experiments. 


Case  A 

Computed 

(air) 

Nominal 

(air) 

Moo 

6.34 

7.39 

Too  (K) 

548.2 

565.3 

Tv oo  (K) 

851.4 

565.3 

Poo  (g/m3) 

47.77 

15.25 

tioo  (m/s) 

2974.7 

3521.7 

Table.  8  Computed  free-stream  conditions  at  the  exit-plane  for  case  A. 

7.  Conclusions 

Detailed  simulations  of  hypersonic  high  enthalpy  experiments  in  the  LENS  expansion 
tube  were  performed.  We  focused  on  one  of  the  two  experimental  conditions  provided. 
We  compared  our  simulation  results  with  measurements  made  during  the  experiment, 
in  particular,  Pitot  pressure  and  wall  pressure  at  the  exit  of  the  expansion  tube.  The 
comparisons  showed  that  the  CFD  predicts  correctly  the  arrival  of  the  test  gas  in  the  nozzle 
behind  the  primary  shock.  We  were  unable  to  obtain  complete  quantitative  agreement 
between  the  measurements  and  simulation  results  for  this  experiment.  The  pressure  near 
the  wall  at  the  end  of  the  expansion  section  was  in  good  agreement,  but  the  Pitot  pressure 
along  the  centerline  of  the  tube  at  that  location  was  not  matched.  This  discrepancy  is 
probably  a  result  of  numerical  dissipation  and  insufficient  grid  resolution. 

We  obtained  an  estimate  of  the  available  test  time  during  the  experiment  using  this 
approach.  Under  the  present  conditions,  the  available  test  time  is  about  3  msec.  We 
computed  conditions  at  the  exit-plane,  and  these  can  be  used  as  free-stream  conditions  for 
the  double-cone  flow.  Also,  the  free-stream  resulted  from  this  experiment  is  non-reacted 
and  exhibits  only  a  small  degree  of  vibrational  nonequilibrium.  The  simulation  showed 
that  the  flow  is  uniform  in  a  relatively  large  region  at  the  center  of  the  nozzle  when  it 
enters  the  test-section.  This  region  is  large  enough  to  ensure  that  there  is  uniform  flow 
around  the  entire  double-cone  model.  However,  under  the  present  conditions,  we  expect 
that  the  double-cone  flowfield  will  likely  be  transitional.  Finally,  we  have  demonstrated 
that  it  is  possible  to  simulate  the  entire  experiment  in  this  large  facility  with  reasonable 
accuracy.  This  process  will  help  design  new  expansion  tube  experiments  with  realistic 
free-streams. 
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Summary  and  Conclusions 


1.  Summary 

This  thesis  discusses  numerical  simulations  of  experiments  of  hypersonic  flows  over 
double-cone  geometries  performed  at  the  Large  Energy  National  Shock  (LENS)  facility. 
This  work  was  motivated  by  earlier  computational  studies  of  double-cone  and  double¬ 
wedge  experiments,  which  showed  primarily  that  flows  with  strong  shock  interactions  may 
serve  as  test  cases  for  the  validation  of  nonequilibrium  chemistry  models.  The  double¬ 
cone  flow  was  chosen  because  it  exhibits  circular  symmetry,  thus  only  two-dimensional 
simulations  are  required,  which  substantially  decreases  the  cost  of  such  simulations. 

First,  the  governing  equations  for  the  relevant  flow  physics  were  presented,  along  with 
the  physical  models  that  provide  closure,  namely  the  transport  coefficients,  the  chemical 
reaction  rate  coefficients,  the  equilibrium  constants,  and  the  vibration-dissociation  coupling 
model.  Also,  the  discretization  of  the  governing  equations  and  the  numerical  method  of 
solution  were  discussed  in  detail. 

In  the  first  stage  of  this  work,  simulations  of  double-cone  flow  experiments  at  low  en¬ 
thalpy  conditions  were  performed,  with  the  purpose  to  assess  the  ability  of  the  numerical 
method  to  reproduce  the  basic  flows.  These  experiments  were  performed  at  low  pressures, 
such  that  the  flow  over  the  model  was  laminar.  There  was  no  prior  knowledge  of  the  ex¬ 
perimental  measurements  when  these  simulations  were  performed.  The  blind  comparisons 
showed  good  agreement  between  the  numerical  predictions  and  the  experimental  surface 
measurements,  with  the  exception  of  the  heat  transfer  rate  at  the  first  cone. 

In  order  to  better  understand  the  limitations  of  the  numerical  method  and  to  assess 
its  accuracy,  the  double-cone  flow  was  simulated  using  different  flux  evaluation  methods  of 
varying  degree  of  sophistication.  In  this  study,  the  effect  of  numerical  dissipation  exhibited 
by  different  flux  evaluation  schemes  was  examined.  Additionally,  grid  convergence  studies 
were  presented.  As  a  result  of  these  studies,  we  established  a  measure  for  the  numerical 
dissipation  and  the  rate  of  grid  convergence  of  a  given  numerical  scheme,  based  on  the 
computed  size  of  the  separation  zone  that  is  resulted  by  employing  this  scheme. 

Next,  the  reason  for  the  discrepancies  observed  in  simulations  of  low  enthalpy  flows 
was  investigated  by  performing  detailed  simulations  of  the  experiments  and  the  wind  tunnel 
flow.  We  used  CFD  to  simulate  the  flow  in  the  nozzle  of  the  LENS  tunnel.  By  matching 
Pitot  pressure  measurements  of  survey  instruments,  we  attempted  to  compute  the  free- 
stream  conditions  at  the  exit-plane  of  the  nozzle.  We  found  that  there  is  some  degree 
of  nonequilibrium  in  the  free-stream  due  to  the  low  density  at  which  the  experiments 
were  performed.  This,  along  with  weak  vibrational  energy  accommodation  at  the  surface 
is  responsible  for  the  discrepancies.  From  this  approach,  a  methodology  was  developed 
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for  designing  and  characterizing  the  LENS  experiments.  Following  this  study,  new  high 
enthalpy  experiments  were  performed.  We  applied  this  methodology  to  the  new  high 
enthalpy  double-cone  flows  of  nitrogen  and  air.  We  observed  good  agreement  for  nitrogen 
flows,  but  we  were  unable  to  reproduce  the  measurements  for  air  flows. 

Finally,  an  alternative  method  for  producing  high  enthalpy  non-reacted  ffee-streams 
was  presented.  This  makes  use  of  the  newly  developed  expansion  tunnel  facility  at  LENS. 
With  this  study,  we  showed  that  useful  test  times  of  up  to  3  msec  can  be  obtained,  and 
the  free-stream  exhibits  only  a  small  degree  of  vibrational  nonequilibrium,  but  maintains 
its  original  composition. 

2.  Conclusions 

This  work  showed  that  computational  fluid  dynamics  is  capable  of  accurately  repro¬ 
ducing  hypersonic  flows  with  shock/shock  and  shock  wave/boundary  layer  interactions. 
However,  large  computational  meshes  may  be  required,  and  it  is  important  to  perform  a 
careful  grid  convergence  study.  Also,  for  these  flows  with  large  regions  of  separation,  the 
amount  of  physical  time  that  must  be  simulated  may  exceed  several  flow  times.  Thus, 
obtaining  accurate  results  on  a  fine  mesh  requires  the  use  of  an  implicit  solver. 

The  double-cone  flow  is  an  excellent  case  for  validating  computational  fluid  dynamics 
codes.  This  is  in  part  because  the  embedded  region  of  separation  is  sensitive  to  high 
enthalpy  effects,  and  its  size  can  be  measured  accurately  during  the  experiment  from  the 
surface  measurements.  The  LENS  experiments  use  densely  instrumented  models,  and  allow 
for  an  accurate  measurement  of  the  separation  zone  size  to  be  made.  Thus,  provided  that 
a  CFD  code  accurately  predicts  the  flow  under  conditions  where  chemical  reactions  are 
not  taking  place,  air  chemistry  models  can  be  tested  at  higher  enthalpy  conditions. 

Also,  double-cone  flow  predictions  are  sensitive  to  the  quality  of  the  numerical  method 
employed.  In  particular,  the  predicted  size  of  the  separation  zone  on  a  given  computational 
mesh,  is  directly  related  to  the  amount  of  numerical  dissipation  in  the  scheme.  In  general, 
less  dissipative  numerical  schemes  predict  a  larger  separation  zone  than  the  more  dissipa¬ 
tive  ones.  Additionally,  dissipative  numerical  methods  may  require  many  more  points  to 
produce  grid  resolved  results.  Therefore,  we  can  quantify  the  amount  of  numerical  dissi¬ 
pation  in  a  given  scheme  by  the  separation  zone  size  that  is  predicted  using  the  scheme. 

In  validation  studies,  it  is  crucial  to  have  the  correct  free-stream  specification  and 
boundary  conditions,  but  also  adequate  modeling  of  the  flow  physics.  For  the  low  density 
flows  of  interest,  it  is  important  to  model  the  gas  surface  interaction  correctly,  and  account 
for  non-continuum  effects.  In  doing  this,  partial  accommodation  to  the  wall  of  the  internal 
energy  modes  of  the  gas  must  be  considered. 


96 


Because  the  separation  zone  size  depends  on  the  grid  resolution  and  the  quality  of  the 
numerics,  it  is  possible  to  obtain  spurious  “agreement”  with  experiments  due  to  poor  grid 
resolution.  Therefore,  the  double-cone  simulations  illustrate  the  importance  of  conducting 
careful  grid  convergence  studies.  Also,  for  the  same  reason,  we  see  the  value  in  making 
blind  comparisons  with  experiments. 

3.  Recommendations  for  Future  Research 

The  development  of  double-cone  experiments  at  high  enthalpy  remains  a  difficult  task. 
The  process  of  specifying  the  experimental  conditions  remains  a  source  of  uncertainty,  par¬ 
ticularly  for  air  flows.  This  is  because  the  chemical  reaction  rates  are  not  accurate  for  the 
sudden  expansion  of  high  temperature  air  that  is  partially  dissociated  and  in  nonequilib¬ 
rium.  Thus,  we  do  not  have  complete  knowledge  of  the  state  of  the  gas  in  the  hypersonic 
nozzle  and  at  the  exit-plane.  However,  using  CFD  to  understand  the  nozzle  flows  has 
provided  an  advantage  over  traditional  methods,  and  more  work  needs  to  be  done  in  that 
area.  This  process  must  be  aided  by  experimental  measurements  made  in  the  nozzle. 

Using  the  expansion  tube  for  generating  hypersonic  free-streams  at  high  enthalpy  can 
aid  the  validation  process.  First,  experiments  must  be  designed  to  produce  low  density, 
low  enthalpy  free-streams,  resulting  in  essentially  perfect  gas  laminar  flows.  In  this  way, 
comparisons  across  facilities  can  be  made.  Upon  successful  completion  of  this  process,  high 
enthalpy  experiments  should  be  designed  and  characterized  with  CFD. 

Finally,  independent  work  on  high  temperature  nonequilibrium  air  chemistry  models 
should  be  done.  This  is  an  area  of  research  where  work  in  other  disciplines,  such  as 
computational  chemistry,  must  be  done. 
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